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Section  I 


INTRODUCTION 

The  determination  of  the  skin  currents  induced  on  a mis- 
sile in  flight  is  an  important  first  step  in  predicting  the 
response  of  internal  subsystems  to  an  impressed  electromagnetic 
field.  The  missile  skin  currents,  however,  may  be  affected 
by  the  presence  of  the  rocket  exhaust  (plume)  which  consists 
of  partially  ionized  gases  or  plasma  and  therefore  has  a 
finite  conductivity.  Because  of  the  difficulty  and  expense 
involved  in  obtaining  measured  data,  a reasonably  accurate 
numerical  solution  procedure  is  desirable.  Furthermore,  num- 
erical results  can  provide  a partial  validation  of  actual 
experiments.  The  more  sophisticated  numerical  analyses  gen- 
erally involve  a body  of  revolution  rcodel  for  both  the  mis- 
sile and  the  plume.  In  this  report  we  consider  body  of 
revolution  models  in  general,  for  both  perfect  electric  con- 
ductors and  dielectric  bodies,  without  reference  to  a parti- 
cular structure. 

The  primary  impetus  for  the  development  of  the  computer 
code  described  in  this  work  has  been  the  desire  for  simplicity 
of  formulation  and  for  increased  accuracy  and  efficiency. 

Other  numerical  formulations  for  bodies  of  revolution  are 
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already  available  [1,  2,  3]  and,  indeed,  the  formulation  of 
[3]  has  been  generalized  to  the  missile/plume  problem  and 
results  have  been  presented  in  a recent  report  (Wu,  et . al. 
[3]).  However,  some  stability  problems  have  been  encountered 
in  these  formulations  for  certain  geometries,  whereas  the 
formulation  presented  in  Section  II  has  been  extensively 
tested  and  has  demonstrated  no  apparent  stability  problems. 
More  importantly,  however,  this  formulation  and  computer  code 
serves  as  a basic  building  block  from  which  an  extremely 
complicated  code  is  to  be  developed  which  treats  a very  sophis 
ticated  model  -,r  the  missile/piume  structure  in  which  the 
inhomogeneous  plume  is  modeled  by  layers  of  homogeneous  mat- 
erial. This  formulation  will  be  described  in  a forthcoming 
report . 

The  numerical  solution  procedure  for  a body  of  revolution 
which  may  be  either  a perfectly  conducting  body  or  a dielectri 
body,  is  presented  in  Section  II.  Numerical  results  for  sev- 
eral cases  are  also  presented  ar.d  discussed.  In  Appendix  B 
the  implementation  of  the  numerical  formulation  is  discussed 
with  direct  reference  to  an  available  computer  code  which  if 
also  listed  In  the  appendix. 
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Section  II 


NUMERICAL  SOLUTION  PROCEDURE 
FOR  A BODY  OF  REVOLUTION 

In  this  section  integral  equations  are  derived  for 
equivalent  surface  currents  induced  on  a dielectric  body 
of  revolution  subject  to  plane  wave  illumination.  Scatter- 
ing by  a perfectly  conducting  body  of  revolution  can  also 
be  obtained  as  a special  case.  The  numerical  solution  pro- 
cedures are  described  and  the  method  of  moments  is  applied. 
Numerical  results  are  obtained  for  several  cases  and  are 
compared  with  other  available  data. 

2 . 1 Formulation  of  the  Integral  Equations 

Consider  the  body  of  revolution  of  Fig.  2.1,  where 
the  body  is  formed  by  rotating  a planar  curve,  called  the 
"generating  arc,"  around  the  z-axis  (also  called  the  axis 
of  the  body  of  revolution).  The  regions  exterior  and  in- 
terior to  the  body  are  denoted  as  regions  1 and  2,  respec- 
tively. The  t-coordinate , which  is  depicted  in  Fig. 2.1, 
follows  the  generating  arc  on  the  body  surface  S.  The  body, 
with  constitutive  parameters  (p2>  e2,  » 0) , is  considered 
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Figure  2.1.  Geometry  of  the  body  of  revolution. 


to  be  Immersed  in  an  infinite  homogeneous  medium  with  pa- 
rameters (Pp  Cp  o^O).  The  generalization  to  non-zero 
conductivity  of  either  the  body  or  exterior  medium  param- 
eters is  trivial,  but  has  not  been  considered  here. 

Boundary  conditions  require  that  the  total  fields  tan- 
gential to  the  surface  of  the  body  be  continuous  from  Region 
1 to  Region  2.  These  conditions  may  be  expressed  as 


ft  x (E  - E8)  - 

ft  x 

EinC 

(2 

. la) 

ft  x (H  - H8)  - 

•\ 

ft  x 

HlnC 

(2 

.lb) 

where  (Einc,  Hinc)  constitute 

the 

incident  fields, 

<E8, 

H8) 

are  the  scattered  fields  in  Region  1,  and  (E,  H)  are  the 
total  fields  in  region  2,  and  n ■ $ * t is  the  outward  unit 
normal  to  the  body  surface.  Via  the  equivalence  principle, 
the  body  is  replaced  by  two  sets  of  electric  currents  ( J ^ , 
i “ 1,2)  and  magnetic  currents  (M^,  i»l,2)  — one  set  just 
inside  the  surface,  and  the  other  just  outside  the  surface. 
Each  of  ..vs  two  sets  of  currents  radiates  in  an  infinite 
homogeneous  medium  having  the  constitutive  parameters  with 
the  correspoivding  medium  index,  i ■ 1,2.  Thus  the  fields 
Indicated  in  (2.1)  may  be  expressed  as 

Es(r)  - -jwA1  (r)  - Vt^r)  - --  V x Fj  (r)  (2.2a) 


5 


(2.2b) 


HS(r) 


■jwFj,  (r) 


E(r)  - -juA2(r)  - 


L (r)  + “ V * At(r) 


V$2 ( r ) - “ V * F2(r) 


(2.2c) 


H(r)  - -jioF2(r) 


- Vf2(r)  + “ V * A2(r)  , 


(2. 2d) 


where  the  potentials  are  defined  by 


Vr)  ■ vF 


If 


Ji(r,)Gi(t,  r')  dS' 


(2.3a) 


J/5 


Fi(r)  “ 4i  I |Ml(7,)Gi(7’  r'>  dS' 


(2.3b) 


>i<')  ' 4^7  r’ 


) dS’ 


(2.3c) 


Vr)  ■ x™  1 1 pi<r,>Gi(r'  r,)  ds’  * (2,3d) 


i ■ i *2  , 


and  where 
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G1(r,  r’) 


e 


(2.4a) 


~R  » 

i - 1,2  , 

R * | r — r * J » j^p ^ + p*^  — 2pp*cos(4>  — 4>*)  H-  (z  — 3s*)^| 

(2.4b) 

The  continuity  of  the  tangential  fields  at  the  boundary 
requires  that 

*1  m E 1 * (2.5a) 

HjL  - -M2  S M . (2.5b) 

Note  that  we  have  introduced  unsubscr ipted  currents  J and 
M.  One  may  then  also  define  an  unsubscripted  charge  by 

p*  = p*  * -p*  - J[VJ«J(r')]  (2.6a) 

p“  = P*  * -p“  * £[VJ*M(r')]  . (2.6b) 

Combining  Eqs.  (2.1)  through  (2.6)  allows  one  to  write  two 
vector  integro-dif f erential  equations  which  may  be  enforced 
on  the  surface  of  the  body  to  obtain  the  unknown  electric 
and  magnetic  currents: 
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a x iinc 


a x ' 4¥ 

47T 


J (y1G1  + p 


2G2>  dS’ 


+*v  //(v«,j)[^+y ds* 


//- 


+ 47  7 * I I M (G1  + G2)  dS’ 


n x HlnC 


a 

9 *{ ^ JJ5  <E‘G> 


+ e2G2)  dS' 


//":-»Et  * g 


dS 1 


(2.7a) 


// 


4r  7 « II  J (Ct  + G2)  dS' 


(2.7b) 


where  the  dependence  of  the  appropriate  quantities  on  source 
and/or  field  coordinates  is  understood.  Since  the  term  in- 
volving the  curl  operator  appearing  in  (2.7)  is  not  contin- 
uous at  the  boundary,  the  direct  interchange  of  integration 
and  differentiation  is  not  allowed.  It  can  be  shown  [4], 
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however,  that  on  the  surface  S, 


V x 


// 


U (G j|  + G2)  dS'  - - 


ti 


U X V(G1  + G2)  dS' 


, (2,8) 

U « J or  M , 


where  j-f  represents  a Cauchy  Principal  Value  integral. 

S 

Eqs.  (2.7)  can  be  written  in  a component  operator  form  as 


E‘"C  • BU<V  +B,j(J+)  +B„(Mt)  +Bu(M,)  (2.9a) 

EJ"C  ' 621<V+  622<V+e23(Mt)  + 824(V  <2'9b> 

"J”'  ' 831<Jt)+  832'V+  633(Mt>  + 834(V  <2-9<=> 

' 841<Jt)+  842(V + 643(Mt)  + 844<V  ’ (2'9d) 


where  is  the  appropriate  integro-dif f erential  operator. 
It  is  desirable  to  express  all  of  the  quantities  of  (2.9) 
in  terms  of  the  local  coordinates  (t,<P)  on  the  body  surface. 
An  orthogonal  triad  of  unit  vectors  (n,$,t)  may  also  be 
associated  with  each  coordinate  point  ( t , 4> ) , where  n, 

A 

and  t are  defined  as  follows: 

n •*  cosycos<}>  x + cosY8in4>  y - siny  z (2.10a) 

AAA 

<p  « -sin<{>  x + co8^  y (2.10b) 

t » sinYco8<{»  x + sinYsimfi  y + cosy  z , (2.10c) 
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where  y is  the  angle  between  the  tangent  to  the  generating 

A /\ 

curve,  t,  and  the  z-axis,  defined  to  be  positive  if  t points 

A 

away  from  the  z-axis  and  negative  if  t points  toward  the 
z-axis.  We  note  that  the  surface  divergence  in  this  co- 
ordinate system  becomes 


V'  *U 
s 


_1_ 

P* 


tr<V 


(2.11) 


U - J or  M . 


One  may  now  expand  (2.7)  into  components  and  compare  to 
(2.9)  to  obtain  expressions  for  the  : 


6n<V 


.1“ 

4tt 


SS 


Jfc  [sinYainy ' cos(4>  ’ -$)+cosyco8y  * ] { y ^ G j -4-y 2 G 2 


ds 


+ 


_X_  i_ 

4noj  3t 


(2.12a) 


g12(V 


4tt 


// 


J^sinysin^'-^Hv^G^+y^G^  dS' 


4tto) 


9_ 

3t 


Si 


i a 


w (V 


> dS  ' 


(2.12b) 
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8,3<Mt> 


[ (p  ' sinycosY ' - pcosYsinY^sin^’- 


S 


+ (z-z  ' ) sinYflinY  ' ain((}i ' -<j))  ] ^{G^+G2)  ' 


®U<V 


’cosy  - pcosYcos  (<p  ' -$) 


S 


+ (z-z  ' ) sinYcos  ($  ' -<j>)  ] ^-{G^+Gj,}  ^S' 


621(Jt)  "4? 


If 


Jt8inY,sin((})’-4'){M1G^+p,G2}  dS’ 


622(>V 


dS  ’ 


+ 


S 


!r(V 


dS’ 


(2.12c) 


(2 . 12d) 


( 2 . 12e) 


(2. 12f ) 
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823 (Mt) 


4ir  JJ  R 


^[pcosy*  - p ' cosy ' cos  (<j) ' -<{>) 


(z-z  ' ) 8iny  ’ cos  (4> ' -<l>)  ] ^r^gi+G2^  ds* 


(2 . 12g) 


B24(V- 


JL_ 

4ff 


ff 


M,  ,, 

-jH~(z~z  * )ain(<}>  ’ -$)  — {Gj+G^  ds' 


(2 . 12,h) 


where  R is  defined  by  (2.4)i  and 


dGi  (l  + jkiR)  -jktR 

___  - e 

i - 1,2  • 


(2.13) 


Eqs.  (2.12)  define  the  3^  for  (2.9a)  and  (2.9b)  . Recall 
that  3^  (U)  is  an  operator  on  the  function  U.  One  may  also 
consider  8^,  i*l»2,  J — 1,2,  to  be  dependent  on  the  param- 
eters of  the  media,  i.e. 

B i j (u ) “ Bi:J  (U ; Uj,  Ej,  u2,  e2)  , 

i - 1,2  , 

J - 1,2  . 

Similarly , 

8^j(0)  * 8 i ^ (U ; p ^ ^ » y 2 » e2^  * 

i - 3,4  , 
j - 3,4  . 
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One  then  finds  that 


33(Mt; 

V 

V 

y2. 

e2> 

■ 8U<V 

V 

V 

e2> 

u2) 

(2.14a) 

34  (M4> ; 

V 

ei* 

y2. 

£2> 

' 812<V 

V 

V 

e2  * 

V 

(2.14b) 

43(Mt; 

V 

ei> 

V 

e2> 

' 821<V 

V 

e2  * 

V 

(2.14c) 

44 <V 

V 

V 

U2‘ 

e2} 

' S22(V 

\ 

ei» 

V 

e2» 

u2) 

(2 . 14d) 

and 


* -813(Jt> 

( 2 . 14e) 

32 ( 

* -8h<V 

(2. 14f) 

' '823(Jt) 

( 2 . 14g) 

M<V 

" "824 

(2 . 14h) 

Thus  Eqs.  (2.14),  which  are  actually  statements  of  duality 
[ 5 ],  serve  to  define  the  in  (2.9c)  and  (2.9d)  via 
Eqs.  (2,12).  Furthermore,  if  the  body  is  a perfect  conduc- 
tor, (2,9)  reduces  to 


inc 

4> 


0ll(Jt;Wl,el,y2“°’e2“OO)  + ei2^J4>;Ml,el,U2"0,e2*“) 

(2.15a) 

+ 822(VW‘,2'0-e2”’)  ' 

(2.15b) 
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For  a numerical  solution  of  (2.9)  or  (2.15) the  gener- 
ating arc  is  approximated  as  a sequence  of  .Linear  segments 
as  depicted  in  Fig. (2. 2)  where  the  approximation  to  the 
generating  arc  Is  shown  in  the  plane  <J>  « 0.  This  segmented 
generating  arc  is  rotated  about  the  z-axis  to  obtain  an 
approximation  to  the  surface  of  the  body  of  revolution. 

The  points  t^,  t ^ , ....  tN+1  specify  the  end  points  of  the 
linear  segments  approximating  the  generating  arc  and  are 
written  in  terms  of  the  coordinates  p and  z.  The  "half 
points"  t^,  t^»  •••»  are  defined  as 

(t  + t . ) 

V*  ■ ' 2 - (2-16) 

1 < n < N+l  . 


The  variations  of  the  unknown  electric  and  magnetic  currents 
flowing  on  the  surface  are  approximated  by  pulse  functions 
in  the  t-direction  and  are  expanded  in  Fourier  series  in 
the  (^-direction.  The  expansion  of  the  electric  current  is 
given  by 


j(t',<r) 


» N+l 


+ i >2 

m*-60  n*l 


(2.17a) 


] 4 


Figure  2.2.  Approxiaation  of  the  generating  arc  by  linear  segments. 
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The  derivative  of  with  respect  to  t,  which  contributes 
to  the  charge,  is  approximated  as 


m*-“  n*l 


(2.17b) 


In  the  preceding, 


P“(f  ) - 


p”(f) 


1 • V*  S C'  s tn+>i 
0 , otherwise 


[l  , t , < t'  < t 
J n-L  “ n 

1 0 , otherwise 


(2.18a) 


(2.18b) 


At  * t - t . 
n 1 n n-1 


Up  - p ,)2  + («  - 2 ,)21  * (2.19) 

^ n Kn-1  n n-1  1 ’ 


and  I®n  is  the  coefficient  of  the  "total"  current  flowing 
in  the  t-direction  as  defined  by 


X t ( t » , 4> f ) - 2ttp  ' Jt  (t  * , <J>  ’ ) . (2.20) 

This  "total"  current  has  no  physical  meaning  (except  for 
m»0;  for  mj*0,  the  Integrated  current  vanishes)  but,  as  a 
mathematical  concept,  its  use  is  advantageous  in  the 
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solution  procedure.  In  (2.17b)  it  is  assumed  that 

m,0  _ Tm,N+I  _ n • 

1t  1t  ~ 0 

The  magnetic  current  M(t*  and  its  derivative, 

|tt[p  (t  ’ ' ) ] , are  defined  similarly  with  K®n«2rrp'M™n 

O t t t t 

and  M®n  replacing  I°n  and  J™n,  respectively. 

The  current  expansion  scheme  (2.17)  has  a number  of 
advantages.  The  first,  and  most  obvious,  is  that  the 
Fourier  components  of  the  current  can  be  decoupled;  thus 
one  can  solve  for  each  unknown  Fourier  component  of  the 
current  distribution  independently.  The  remaining  advan- 
tages arise  from  the  use  of  the  staggered  pulse  basis  sets 
indicated  by  (2.18),  To  further  illustrate  some  of  these 
advantages  we  consider  the  case  cf  a perfectly  conducting 
cylinder  (Fig.  2.3)  which  has  one  closed  end  and  one  open 
end.  Note  that  the  "total"  current  flowing  in  the  t-direc- 
tion  (It)  is  zero  for  all  Fourier  components  at  t « 0 and 
t*b  [ 6,7  ] and  should  therefore  be  well  represented  in 
the  neighborhood  of  these  two  points  by  half  pulses  vlth  a 
zero  coefficient,  as  in  the  cases  of  the  TE  strip  and  the 
rectangular  bent  plate  described  in  other  works  [8,9].  On 
the  other  hand,  the  current  density  in  the  ({(-direction  (J^) 
approaches  zero  (u  i1  ±1)  or  a constant  (m*±l)  as  t 0 [ 6,7  ] 
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Figure  2,3.  Cross  section  of  a perfectly  conducting  cylinder  with 
one  open  end. 


18 


i 


and  is  always  singular  at  t*b,  as  required  by  the  so-called 
edge  condition  [10].  It  has  been  demonstrated  in  the  pre- 
vious works  [8,9] that  the  staggered  subdomain  scheme,  such  as 
(2.17)  and  (2.18),  employing  full  pulses  to  represent 
adjacent  to  the  points  t*0  and  t*b,  can  accurately  model 
the  singular  currents  near  edges  and  should  clearly  model 
a finite  current  as  well.  The  full  pulse,  however,  may  not 
model  a current  which  approaches  zero  well,  but  in  view  of 
the  fact  that  such  a situation  occurs  only  at  points  where 
the  surface  intersects  the  axis  of  the  body  of  revolution 
(p  • 0)  and  that  the  current  moment  represented  by  the  pulse 
is  therefore  relatively  small,  this  deficiency  should  not 
significantly  affect  integral  functionals  of  the  current 
such  as  radar  cross  section,  etc.  Thus  the  staggered  pulae- 
Fourier  series  basis  set  enjoys  the  same  advantages  de- 
scribed for  the  staggered  subdomain  scheme  on  the  rectangu- 
lar beut  plate,  namely,  the  current  expansion  ensures  that 
t-dlrected  components  of  current  vanish  at  knife-like  edges 
and  are  continuous  at  structure  edges  (sharp  bends),  and 
that  the  ^-directed  current  and  the  charge,  both  of  which 
are  singular  at  an  edge,  are  allowed  to  vary  independently 
on  opposite  sides  of  an  edge.  The  basis  set  (2.17)  there- 
fore allows  one  to  model  both  open  and  closed  bodies,  and 
bodies  with  sharp  edges,  with  no  special  procedures  required 
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at  edges  or  points  where  the  surface  intersects  the  body 
axis.  Of  course,  a dielectric  body  has  no  knife-like  edges, 
but  it  may  have  structure  edges  at  which,  depending  on  the 
parameters  of  the  medium,  it  may  be  necessary  to  represent 
currents  which  are  singular  [11].  The  basis  set  (2.17), 
which  is  used  for  both  the  electric  and  magnetic  currents, 
provides  accurate  modeling  in  the  vicinity  of  such  edges. 

We  next  define  the  testing  functions 

T”<t,#>  - Pj(t)e“jp<,)  (2.21a) 

TPq(t,<j>)  - Pq(t)  e~^P^  . (2.21b) 

Eqs . (2.9a)  and  (2.9c)  are  tested  with  (2.21a),  while  (2.9b) 

and  (2 . 9d)  are  tested  with  (2.21b)  . The  testing  procedure 
results  in  two  surface  integrations  being  required  for  each 
3^. . However,  we  note  that  the  current  has  already  been 
expanded  as  a Fourier  series  in  the  variable  <j>*.  Likewise 
the  various  kernel  terms  appearing  in  (2.12)  may  be  expanded 
in  Fourier  series  of  the  form 


e-jt‘R0  , 


R, 


00 

I 

m»-<n 


°i e 
m 


Jm($  - <f>’) 


i - 1,2  , 


and 


20 
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jm(4>  - (fr*) 


-3ki»0 


00 


Hat— oo 


with  Fourier  coefficients 
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i - 1,2  , 


7T 

L 


-^lR0 


i I R. 

” J-,  0 


co8(mO  d£ 


(2.22a) 
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_1_  d 

R0  dR0 


-^iRol 

k 

Ro 

md 


cos (m£)  d£  , (2.22b) 


where 


Rq  - [p^  + p ’ ^ - 2pp ’ cos( + (s  - z ' ) 2 ] . (2.22c) 


Testing  with  the  testing  functions  (2.21)  permits  all  inte- 
grations in  the  variables  $ and  to  be  carried  out  analyt- 
ically, simultaneously  decoupling  the  equations  with  respect 
to  their  dependence  on  the  Indices  p and  m.  The  double  in- 
tegral in  the  t-direction  is  reduced  to  a single  Integral 
via  appropriate  approximations  (such  as  used  in  [8,9]).  We 
also  note  that  derivatives  with  respect  to  <j>  and  <j>'  can  be 
performed  analytically  during  the  testing  procedure.  By 
means  of  such  approximations  and  procedures,  one  obtains 
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of  the  generalized  imped- 


expressions  for  the  elements  8^ 

J m 

ance  matrix,  where  the  subscript  m refers  to  the  Fourier 
coefficient  index  and  the  superscripts  q and  n refer  to  the 
Indices  of  the  field  point  and  the  current  source,  respec- 
tively. For  convenience  in  writing  these  expressions , we 
first  define  two  frequently  used  integral  functions: 


f2 

'1<ti»t2;vB>  * J Gi  (v*') 

ID 


dt 1 (2.23a) 


t 

I 


G (t  ,t')p'  dt' 
m q 


(2.23b) 


where  is  defined  by  (2.22),  We  also  introduce  the  aux- 
n 

iliary  "weight"  functions  arising  from  testing  in  the  t- 
direction : 


X (At  , y ) 
As  q 'q 

X (At  ,Y  ) 

Ac  q q 


ata+l,11,y<1+l+tto,lny<1 

2 

At  ..cosy  ,,+At  cosy 

ail Li±I____g Li 

2 


(2.24a) 


(2.24b) 


The  expressions  for  the  elements  of  the  generalized  imped- 
ance matrix  are  given  on  the  following  pages: 
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The  U's  in  (2.25)  are  integral  functions  defined  by 


ml  11 

VVW>-t  / g- dtdf 

•'t,  *'-71  0 ° 


(2.26a) 


t2  tr 


Ul^tl,t2;tq’m^ 


H 

•'tj  •'-TT 


~ Z">  sin  (mO  sing  d£dt'  (2.26b) 


C2  W 

Ul(tl,t2;tq*m)  " jf  f p'd^dt' 


U2^tl,t2;tq*m^ 


t„  n 

a. 


(2.26c) 


.^-r.z,)coa(mf-)coag  ijjL  dUt.  (2 . 26d) 


u2(t 


t2  U 

l»t2;tq»m^  * £ ^ r ~~~^c  °s(mg)co8g  p'd^dt' 


U3(tl*t2Stq*m> 


tj  IT 

T I j^"1 8in(mOsin^  dSdt 

J,  J./o  dR0 


(2. 26e) 


(2 . 26f ) 


t2  Tf 

U4^tl,t2’tq,m^  “ £ -^P— ^-.P.  ^cos(mC) 


dG 

dR, 


p'dCdt’  (2 . 26g) 


The  integrals  defined  by  the  U's  are  actually  recombinations 
of  expressions  involving  the  Fourier  coefficients  (2.22b) 
and  are  expressed  in  the  manner  shown  to  facilitate  the  sin- 
gularity analysis  of  self  terms,  which  appears  in  Appendix  A. 
We  have  now  completely  defined  the  elements  of  the  impedance 
matrix  corresponding  to  the  operators  appearing  in  (2.9a) 
and  (2.9b);  the  matrix  operators  arising  from  (2.9c)  and 
(2.9d)are  obtained  from  the  functional  relationships 
(2.14).  One  should  note  that  many  of  the  integrals  which 
must  be  computed  in  the  calculations  of  the  elements  (2.25) 
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appear  in  several  pi ices  and  hence  can  be  used  several  times 
to  increase  efficiency. 

The  incident  plane  wave  fields  are  expressed  as 


— *nc  iAi  iAi 

Einc  ■«  (E^e1  + gJOe 

/s  — 

h1"'  - IcEje1  - E«l).'ik,"’r 
n <p  0 


where 


a i 
01 


n 

r 


1.  1 A •£  ^ A £ A 

cos0  co s<)>  x + co80  sin$  y - sin0  z 

i A.  I 

-sin<|>  x + cosifi  y 

-sinS^cosili1-  x - sinB^sin^1  y - cosB*  z 
pcos<()  x + psin<)>  y + zz 


(2.28a) 

(2.28b) 


(2.29a) 
(2.29b) 
(2.29c) 
(2 . 29d) 


and  n is  the  free  space  impedance.  Mote  that  n is  in  the 
direction  of  propagation.  The  elements  of  the  forcing  func- 
tion vector  are  thus  determined  by  finding  the  components 
tangential  to  the  surface  and  testing  with  (2.21).  With 
reference  to  (2.9) » these  elements  are  gxven  by  the  ex- 
pressions on  the  following  page,  where  J (u)  is  the  Bessel 
function  of  the  first  kind  of  order  m.  With  regard  to  the 
expressions  (2.30),  we  comment  that,  with  no  loss  in  gener- 
ality, one  may  set  <j>*  « 0. 
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(poe*z)  • T+N 


A full  solution  for  the  currents  on  the  body  of  revo- 


lution consists,  in  general,  of  an  Infinite  number  of 
Fourier  components.  In  practice,  however,  the  series  is 
truncated  such  that  -M£m<[M,  where  M is  the  maximum  number 
of  positive  Fourier  components  to  be  calculated.  Further- 
more, it  can  be  shown  that  the  current  solution  for  -m  is 
related  to  the  solution  for  +m,  so  that  it  is  only  necessary 
to  compute  the  impedance  matrix  and  drive  vector  for  m«0,l, 
2,...,M.  This  relationship  is  determined  by  comparing  the 
signs  of  the  impedance  matrix  and  a decomposed  drive  vector 
for  the  positive  and  negative  Fourier  component  indices. 

If  one  has  available  the  impedance  submatrices  3.  for  pos- 
itive  m,  then  the  impedance  matrix  for  the  corresponding 
negative  Fourier  component,  -m,  is  given  by 
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The  submatrices  constituting  the  drive  vector  can  be  re- 
solved as  follows: 
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(2.32) 


where  the  superscripts  6 and  $ refer  to  the  components  of 
the  incident  field  arising  from  Eg  and  E* , respectively. 
We  then  find  that 
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Thus  the  decomposed 

vectors  are  related 

solution  vectors  by 
s 

r 


negative  Fourier  component  solution 
to  the  corresponding  positive  component 


(2.34a) 


(2.34b) 


jhich  may  be  combined  to  obtain 


(2.35a) 
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(2.35b) 


Thus  Eqs.  (2.9)  (or  (2.14)  for  a perfect  conductor) 
constitute  a linear  system  of  equations  for  the  current 
coefficients  in  (2,17).  The  Fourier  components  of  the  cur- 
rent in  the  problem  are  decoupled,  however,  and  one  may 
solve  for  each  Fourier  component  current  distribution  in- 
dividually, using  the  impedance  matrix  elements  given  by 
(2.25)  and  the  elements  of  the  forcing  function  vector, 
(2.30).  The  forcing  function  vector  is  then  decomposed  into 
two  drive  vectors  ae  indicated  by  (2,32)  so  that  one  may  ob- 
tain both  the  positively  and  negatively  Indexed  Fourier  com- 
ponents of  the  current  distributions  simultaneously,  through 
the  relationships ( 2 . 35)  . 


Numerical  Results  for  the  Bodv  of  Revolution 


The  numerical  procedures  described  in  the  previous 
section  have  been  incorporated  into  a computer  code,  herein 
referred  to  as  "DBR,"  and  numerical  results  for  the  surface 


currents  on  several  structures  have  been  obtained  and  com- 


pared with  other  available  data.  The  code  is  capable  of 
treating  both  dielectric  and  perfectly  conducting  bodies; 
in  the  latter  case,  the  bodies  may  be  either  open  or  closed 
structures . 

All  of  the  structures  considered  in  this  section  are 

assumed  to  be  illuminated  by  a plane  wave  incident  along 

the  z-axis  with  an  electric  field  polarized  in  the  x-direc- 

tion.  Propagation  may  be  in  either  the  positive  or  negative 

z-directlon.  For  a field  incident  along  the  axis  of  the 

+ i d> 

body  of  revolution,  only  the  Fourier  components  with  e~jr 
variation  (m-il)  are  excited.  Using  the  relationships 
(2.35)  the  surface  currents  can  then  be  written  in  the  sim- 
ple form 


J(t,$)  - J®(t)[ej<l>  + e_;j<,,3  + J^(t)[ej(j)  - e"J<h 

■ J.  (t)cos<j>  + j J , (t)  Bin$  (2.36a) 

t 

M(t,4>)  « M®  (t)  [ej<<>  - e"^]  + M®  ( c)  [eJ*  + e~j(h 

* jMt  ( t)  sinij)  + ( t ) co8<p  , (2.36b) 

where 

J (t)  - 2 J®  ( t ) 

P P 
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The  surface  currents  without  superscripts  in  (2.36)  there- 
fore represent  the  variation  in  the  t-direction  of  the 
actual  surface  current  density  on  the  body,  consisting  of 
a sum  over  the  excited  Fourier  components  and  observed  in 
either  the  <j>  • 0°  (Jt>M^)  or  the  $*90°  (J^,Mfc)  plane.  The 
figures  in  this  section  illustrate  the  t-variation  of  these 
actual  current  densities  when  viewed  in  the  appropriate 
plane.  We  first  consider  results  obtained  for  dielectric 
bodies . 

In  Figs. 2.4  and  2.5  the  electric  and  magnetic  surface 
current  distributions,  respectively,  are  illustrated  for  a 
dielectric  sphere.  Results  calculated  by  means  of  DBR  are 
compared  with  results  obtained  by  Wu  [12],  The  radius  of 
the  sphere  is  k^a  - 1 and  its  relative  constitutive  param- 
eters are  * 1 and  e^-4.  The  sphere  is  excited  by  a plane 
wave  propagating  in  the  positive  z-dlrection.  In  general, 
the  data  obtained  with  the  computer  code  DBR  agree  well  with 
that  of  Wu.  One  should  note,  however,  that  the  results 
calculated  by  means  of  DBR  have  a slight  "glitch"  near  the 
points  where  the  surface  meets  the  body  axis.  It  is  spec- 
ulated that  the  rapidly  decreasing  radius  in  this  region 
may  necessitate  using  a smoother  current  basis  (e.g.,  lin- 
ear) for  the  total  current  I in  the  vicinity  of  the  points 
where  the  surface  meets  the  body  axis,  where  it  is  presently 
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Wu  (N=20) 
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Figure  2.5.  Magnetic  surface  current  distribution  on  a dielectric  sphere  illuminated  by 
an  axially  incident  plane  wave. 


modeled  by  a zero  magnitude  half  pulse.  This  speculation 
has  not  yet  been  confirmed;  however,  numerous  other  attempts 
to  correct  the  problem  have  been  unsatisfactory.  For  open 
conducting  bodies  where  the  current  density  also  ap- 
proaches zero  near  the  ends  of  the  body,  no  such  "glitch" 
has  been  detected.  The  presence  of  the  "glitch"  should  have 
a negligible  effect  on  field  calculations,  however,  except 
possibly  very  near  the  points  where  the  surface  meets  the 
axis,  since  the  total  current  moment  for  the  pulses  in  this 
region  is  small. 

The  electric  and  magnetic  surface  current  distributions 
for  a dielectric  sphere  with  radius  ,2A  and  * 80  are  shown 
in  Figs.  2.6  and  2.7,  respectively,  and  are  again  compared 
with  results  obtained  by  Wu  [12],  One  observes  in  Fig.  2.6 
that  the  agreement  of  the  results  is  excellent  if  one  ig- 
nores the  (apparently)  non-physical  oscillations  exhibited 
by  the  ({(-component  of  current  in  Wu’s  solution.  Note  that 
the  results  obtained  with  DBR  for  N ■ 30  are  presented  as 
continuous  curves,  rather  than  discrete  points.  We  comment 
that  the  data  points  nearest  the  points  where  the  surface 
meets  the  body  axis  have  been  ignored  and  the  curves  have 
been  extrapolated  in  this  region.  The  data  points  which 
were  ignored,  however,  have  been  plotted  in  the  figure  as 
dots  which  do  not  lie  on  the  continuous  curve.  All  other 
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Figure  2 . 7.  Magnetic  surface  current  distribution  on  a dielectric  sphere  illuminated  by 
an  axially  incident  plane  wave. 


data  points  lie  on  the  curves.  The  results  for  the  magnetic 
current  (Fig.  2.7)  are  also  in  good  agreement.  For  the  <£- 
component  of  magnetic  current,  however,  the  data  obtained 
with  DBR  for  N * 20  exhibit  a slight  downward  bulge  near 
t/  (ttX)  = .075,  whereas  Wu’s  results  show  no  such  behavior. 

The  bulge  disappears  when  more  unknowns  are  used,  as  illus- 
trated by  the  curve  for  N * 30. 

In  Fig.  2.8  numerical  data  are  compared  with  the  exact 
solution  for  the  current  densities  on  a "vacuim  dielectric" 
(e  * 1)  cylinder.  The  exact  solution  is,  of  course,  given 
by 
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The  numerical  data  agree  reasonably  well  with  the  exact  so- 
lution when  one  considers  the  fact  that  the  ^-components  of 
current  must  exhibit  a step-function  jump  in  magnitude  at 
the  edges  of  the  cylinder.  The  ability  of  a computer  code 
to  model  such  a case  accurately  may  be  important  if  one 
wishes  to  consider  objects  of  very  low  dielectric  contrast. 
Figs.  2.9  and  2.10  show  calculated  electric  and  magnetic 
current  distributions,  respectively,  for  the  same  size  cyl- 
inder when  the  dielectric  constant  is  increased  to  e «4. 

r 


We  next  consider  numerical  results  for  two  perfectly 
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Figure  2.8.  Electric  and  magnetic  surface  current  distributions  on  a "vacuum  dielectric 
cylinder  illuminated  by  an  axially  incident  plane  wave. 


Figure  2.9.  Electric  surface  current  distribution  on  a dielectric  cylinder  illuminated 
by  an  axially  incident  plane  wave. 
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conducting  bodies.  Fig.  2.11  illustrates  the  current  dis- 
tribution on  a cone-sphere  structure  illuminated  by  an 
axially  incident  plane  wave  propagating  from  the  cone  tip 
toward  the  spherical  cap.  Numerical  results  obtained  with 
DBR  are  compared  with  results  of  Mautz  and  Harrington  [ 1 ] 
and  Poggio  and  Miller  [13] . The  results  of  Mautz  and 
Harrington  are  obtained  via  an  EFIE  formulation,  while 
Poggio  and  Miller  employed  the  MF1E.  Results  obtained  with 
DBR  are  in  excellent  agreement  with  those  obtained  with  the 
MFIE  and  do  not  exhibit  the  (apparently)  non-physical  oscil- 
lations in  the  current  distributions  obtsined  by  Mautz  and 
Harrington . 

We  also  note  at  this  point  that  there  may  exist  a theo- 
retical relationship  between  the  t-  and  <j>-component s of  cur- 
rent at  the  points  where  the  surface  meets  the  body  axis. 
Such  a relationship  is  of  interest  when  treating  body  of 
revolution  structures.  Consider  the  region  near  the  tip  of 
the  cone-sphere  structure  depicted  in  Fig.  2.11.  The  charge 
density  near  the  tip  for  the  m Fourier  component  is  given 
by 

+!?<v] 

' i[Jt!l(r)+  rTF(Jt)  + 5?<V]  • 
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Mautz  and  Harrington  (N=30)  DBR 

Poggio  and  Miller  (MFIE)  _ (N=27) 


the  cone  tip  toward  the  spherical  cap. 


or 


P 


u>r 


^Jtsin(0/2)  +r|^(Jt)  + 


(2.37) 


where  r is  the  cylindrical  coordinate  variable  representing 
the  distance  from  the  z-axis  and  where,  on  the  conical  sur- 
face, r ■ t sin (0/2)  . The  "total"  charge  is  therefore 


q 


(2.38) 


From  symmetry  considerations,  however,  it  can  be  shown  that 
the  total  charge  at  the  points  where  the  surface  meets  the 
body  axis  is  zero  for  m i 0 . At  this  point  we  rely  on  a 
visual  inspection  of  the  numerical  results  to  infer  that 
for  m ■ 1 1 


lim  r 
t-t-0 


3_ 

at 


<v 


o 


(2.39) 


A rigorous  theoretical  investigation  of  the  fields  in  the 
neighborhood  of  the  cone  tip  is  necessary  to  determine  the 
validity  of  (2.39).  Such  an  investigation,  however,  is 
seriously  hampered  by  the  lack  of  tabulated  zeros  of  the 
associated  Legendre  function,  and  has  therefore  not  been 
pursued  in  this  work.  With  (2.38)  we  then  have  for  m * ±1 


Jtsin(8/2)  - , (2.40) 
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which  relates  the  value  of  the  two  current  components  at 
the  point  where  the  surface  meets  the  body  axis.  The  spe- 
cial case  for  0 = it  and  m=±l  was  developed  previously  by 
Mautz  and  Schuman  [7].  Eq.  (2.40)  may  be  of  some  interest 
for  "conically-tipped"  surfaces  when  o*±l.  Referring  to 
Fig.  2.11  again,  one  should  note  that  results  computed  by 
means  of  DBR  (which  enforces  no  special  conditions  at  the 
body  axis)  do  indeed  satisfy  (within  5%)  the  condition 
(2.40)  both  at  the  cone  tip  and  at  the  point  on  the  spher- 
ical cap  at  the  axis. 

When  the  cone-sphere  structure  is  illuminated  by  an 
axially  incident  plane  wave  propagating  from  the  spherical 
cap  to  the  cone  tip,  the  current  distributions  shown  in  Fig. 
2.12  are  obtained.  Results  computed  with  DBR  are  again  com- 
pared with  those  of  Mautz  and  Harrington  [ 1 ] and  Poggio 
and  Miller  [13].  For  the  ^-component  of  current,  agreement 
with  the  MFIE  approach  is  again  excellent,  except  very  near 
the  cone  tip.  However,  for  the  t-component  of  current,  the 
MFIE  approach  yields  results  which  appear  to  oscillate 
slightly.  Results  obtained  by  means  of  DBR  show  essentially 
no  oscillations.  Also,  the  computed  results  again  satisfy 
the  condition  (2.40)  Doth  at  the  cone  tip  and  at  the  point 
on  the  spherical  cap  at  the  axis. 

Fig.  2.13  illustrates  the  current  distribution  on  a 
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Figure  2.12.  Surface  current  distribution  on  a perfectly  conducting  cone-spnere  s 
illuminated  by  a plane  wave  which  is  axially  incident  and  propagates 
the  spherical  cap  toward  the  cone  tip. 


perfectly  conducting  open-ended  cylinder  subject  to  an 
axially  incident  plane  wave.  The  results  computed  with  DBR 
are  compared  with  those  of  Davis,  which  were  obtained  using 
a set  of  hybrid  integral  equations  [ 14  ] . Davis  apparently 
encountered  stability  problems  in  the  solution  cf  the  usual 
EFIE,  and  developed  the  hybrid  equations  as  a means  to  avoid 
these  problems.  A comparison  of  the  data  in  Fig.  2.13 
should  confirm  that  the  formulation  presented  in  this  work 
apparently  does  not  suffer  from  stability  problems  of  the 
kind  encountered  by  Davis.  This  problem  also  presents  the 
opportunity  to  check  the  a posteriori  correction  for  pulses 
which  represent  singular  currents  near  edges  [15],  The  two 
square  symbols  near  the  edges  of  the  plot  in  Fig.  2.13  rep- 
resent computed  currents  which  have  been  corrected  a pos- 
teriori by  dividing  the  computed  current  magnitude  by  /2, 
These  corrected  values  are  in  excellent  agreement  with  the 
results  obtained  by  Davis,  who  used  spline  functions  con- 
taining the  edge  condition.  On  the  other  hand,  if  one 
wishes  to  calculate  the  fields  using  the  computed  currents, 
it  is  probably  better  not  to  correct  the  value  of  the  edge 
current,  since  the  computed  (uncorrected)  value  should  ac- 
curately represent  the  current  moment  over  the  pulse  region. 


Davis 
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Figure  2.13.  Surface  current  distribution  on  a perfectly  conducting  open-ended  cylinder 
illuminated  by  an  axial ly  incident  plane  wave. 


Section  III 


CONCLUSION 

In  this  report  we  have  presented  a s 
numerical  formulation  for  determining  the 
on  both  perfectly  conducting  and  dielectr 
tion.  The  formulation  has  been  implement 
code  (DBR)  which  is  described  and  listed 
erical  results  have  been  presented  for  th 
several  types  of  bodies  of  revolution  and 
ment  of  the  results  with  other  available 
strated.  It  has  also  been  noted  that  oth 
encountered  stability  problems  with  some 
configurations.  There  has  been  no  eviden 
stability  for  any  structural  geometry  whe 
described  in  this  work  are  employed.  A s 
the  solution  for  the  current  has  been  det 
where  the  surface  meets  the  axis  of  the  b 
However,  this  glitch  is  relatively  small 
significant  effect  on  integral  functional 
such  as  fields,  radar  cross-section,  etc. 
has  been  demonstrated  that  the  techniques 
work  provide  accurate  treatment  of  sharp 


imple  and  efficient 
currents  induced 
ic  bodies  of  revolu- 
ed  in  a computer 
in  Appendix  B.  Num- 
e current  induced  on 
the  excellent  agree- 
data  has  been  demon- 
er  researchers  have 
body  of  revolution 
ce  of  solution  in- 
n the  procedures 
light  "glitch"  in 
ected  at  points 
ody  of  revolution, 
and  should  have  no 
s of  the  current 
Furthermore,  it 
described  in  this 
bends  and/or  knife-like 
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edges  in  the  structure. 

The  numerical  formulation  presented  in  this  report  serves 
as  a basis  for  a very  sophisticated  numerical  model  of  the 
missil'e/p  lume  structure  in  which  the  inhomogeneous  plume  is 
modeled  by  layers  of  homogeneous  material.  This  sophisticated 
numerical  model  now  apperars  to  have  been  successfully  imple- 
mented and  will  be  the  subject  of  a forthcoming  report. 
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, APPENDIX  A 

! 

, SINGULARITY  ANALYSES  FOR  SELF  TERMS 

OF  THE  BODY  OF  REVOLUTION  FORMULATION 

When  calculating  elements  of  the  impedance  matrix  for 
the  body  of  revolution  via  (2.25),  one  should  be  aware  that 
the  integrands  of  the  integral  functions  and  U may  possess 
a singularity  when  the  field  point  is  within  the  source  re- 
gion. In  this  appendix  each  integrand  is  investigated  to 
determine  if  a singularity  is  present  and,  if  so,  a numer- 
ical treatment  of  the  integral  function  is  presented.  In  . 
all  of  the  investigations  we  employ  the  following  coordinate 
parameterization  valid  for  t^^StSt^: 

z . .J.J  + Hcos^  ( A. la) 

P - + £sin  , (A.  lb) 

0 < l < Atj  , 

* where  l ■ t-t^_^.  For  all  self  terms  we  then  have  that 

(z-z’)  « (Jl  - £’)cos  Yj  (A. 2a) 

(p-p')  * (£-S,r)siny^  . (A. 2b) 
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A.  1 Singularity  Analysis  of  the  and  ipP  Integral  Functions 

For  a self  term,  the  ijj  integral  function  of  (2,23a)  may 
be  written 


jf  / 


-JkiRo 


cos(mO  d fi, ? , (A. 3a) 

i - 1,2  , 


where 


R0  - £ (P  - P’)2  + 2pp'  (1  - cosD  + (z  - z’)*]  . (A. 3b) 


» \ 2 ■ 


As  t t ' and  £ -*- 0 , Rq  -*■  0 and  the  integrand  of  (A. 3a)  Is 
clearly  singular.  Thus  we  express  ip  as 


*i 


« I 


(A. 4) 


where 


I 


cos(i() 


(A. 5a) 


I 


2 


i 


~ dCd£' 
R0 


(A. 5b) 
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and  where  the  integrand  of 
Furthermore,  w§  have  that 


I.  is  no  longer  singular, 
i 


* [(p  - p')2  + (z  - z’)2  + 4pp,sin2(£/2)]'s 
- Rx  [1  + g28in2(?/2)]]s 

where 

Rj  * [<P  - P1)2  + (z  - z’)2] 


Thus,  with  a simple  change  of  variables,  I-  may  be  written 

i 

as  [16] 


7T 

2 


I,  « 4 


// 


Rjl  + ejsin2?]*5 


dCdi' 


x, 

J 


RjU  + B2]5* 


K 


^U  + 82]35 


dr 


X 

/ 


^k(32)  dr 


(A. 6) 
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where  K(82)  is  the  complete  elliptic  integral  of  the  first 
kind  defined  by 


TT 


1 

[1  - u2sin2(j)3Js 


d<(> 


(A. 7a) 


and 


R2  = [(p  + p')2  + (z  - z * ) 2 3 2 

6*.  ll££l 

p2  R2 


(A. 7b) 


(A. 7c) 


The  integrand  of  (A. 6)  is  still  singular,  however,  and  near 
the  singularity  at  t*t'  varies  as 


r2 

■’  > 1 

[ J,n  A + Hn(R2)  - 

&n(R1)]  . 

t + f 2P 

(A.  8) 

Only  the  last 

term  is  singular,  so  we  add 

and  subtract 

the 

singular  term 

from  <A.6) 

to  obtain 

12 

T a + Tb 

12J  + *2  , 

(A. 9) 

i 

i i 

v/here 

f 2 

V 

/ ft 
*1 

K(8Z)  + 

dA' 

(A. 10a) 
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(A. 10b) 


/ 


AnO^)  di' 


The  Integral  I.  no  longer  has  a singular  integrand  and  the 
b 1 

integral  I-  can  be  evaluated  analytically  by  using  the 

l± 

parameterization  (A. 1)  as  follows: 


j 


S-niR^  dZ' 


At 

1 


An| £ - £ * | di 


~[  (&2  - Aj)  - (£2  - Z)Zn(Z2  - Z) 


- (£  - Z^ZnCZ  - ip] 


(A. 11) 


The  integrals  I.  and  I,,.  can  thus  be  integrated  numerically, 

b i i 

while  I„  is  given  by  (A. 11),  and  we  have 
i 


♦i  - I,  + + i! 

i i i 


(A. 12) 


We  comment  here  that  for  non-self  terms  it  is  more 
t 
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convenient  to  calculate  as 


♦ i " 1i  + 12 
1 ii 


(A. 13) 


where  I„  is  defined  by  (A. 6),  than  to  compute  ^ directly 

l i 

from  (A.  3).  Following  the  procedures  described  above  it 
is  trivial  to  show  that  the  self  term  of  ^ is  calculated 
via 


* I?  + PI?  + Plo 

i i i 


(A. 14) 


where 


cos  (m£)  - 


> d^dl'  . (A. 15) 


A. 2 Singularity  Analysis  of  the  Integral  Function  U. 


The  definition  of  the  integral  function  Uq  is  given  by 


(2.26a)  as 

i 


2 TT 

"o  •/  f [<i  + Jk,R0i  e'3k‘Eo 

"'A  •'-TT 


+ ( 1 + jk2RQ)  e 


!R°  ] 


d£dA' 


, (A. 16) 
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where  is  defined  by  (A. 3b).  As  stated  in  Section  II,  all 
of  the  integral  functions  U may  be  Cauchy  Principal  Value 
integrals.  Thus  one  may  evaluate  the  integrals  by  direct 
numerical  integration,  if  the  quadrature  points  are  chosen 
in  accordance  with  the  definition  of  the  Cauchy  Principal 
Value  Integral.  On  the  other  hand,  evaluation  of  the  in- 
tegrals in  this  manner  may  result  in  the  subtraction  of  very 
large  numerical  values  of  the  integrands  and  hence  lead  to 
possible  loss  of  precision  in  the  calculations.  To  avoid 
such  problems  we  analyze  the  U integral  functions  in  the 
manner  presented  in  Section  a .1.  Near  the  singular  point 
the  integrand  of  (A. 16)  can  be  shown  to  vary  as 

-2mg2 

f 2 2 2l(3/2) 

[<*-*')  + pVJ 

Therefore  we  write 

U0  - I3  + IA  , (A. 17) 

where 
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I.  = 


2 TT  p 

f ( J slnimUslni  r(l+jkiRo)  e'3*1*0 

•'-IT  L R0 


+ (1  + jk2R0)  e 


"jk 


:Ro]_ 


-ImV 


[(*  - V)2  + pV] 


(3/2) 


2 IT 

If 

*0  * . -rr 


-2m£‘ 


J_, r [(«•  - Jl')2  + pV] 


(3/2)  d^dS/' 


>d£<H' 

(A. 18a) 

(A. 18b) 


The  Integrand  of  I ^ is  non-singular  and  can  be  evaluated 
numerically.  The  integral  1^  can  be  integrated  analytically 
and  is  given  by 


I - < 

4 3 1 

P 

( i - ^)Jtn 

Lpir+  1 

n 

1 (Z  - Z:)2  +pV  ] 

"j 

+ 

a2  - Din 

pu  + 

](Z2  - Jl)2  + p2TT2  J 

- - U2  - l)lnW2  - Z)  > . (A. 19) 


Eq . (A. 17)  then  provides  the  desired  result. 
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A. 3 Singularity  Analysis  of  the  Integral  Functions 
and 

The  integrand  of  the  integral  function  U^,  as  defined 
by  (2.26b),  is  identical  with  the  integrand  of  the  integral 
function  Uq  except  for  the  product  factor  (z-z').  It 
should  therefore  be  clear  from  the  analysis  of  Uq  that  as 
Rq  ■*  0 the  integrands  cf  and  are  proportional  to 

-2m£^ (&  - Jl ' ) cosy  * 

P 2 2 2*1  ^ ^ ^ 

which  is  non-singular  and  hence  can  be  integrated 
numerically . 

A. 4 Singularity  Analysis  of  the  Integral  Functions  U£ 
and 

By  comparison  with  the  integrand  of  U^,  one  can  imme- 
diately see  that  as  Rq  0 the  integrand  of  Uj  approaches 

-2U  - & 1 ) cosy ' 

r 2 2 2*1  ^ ^ ^ ) 

[a  - *•>  + PVJ 

This  function  is  singular  and  we  therefore  attempt  to  com- 
pute U2  as 
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U2  = *5  + h 


(A. 20) 


where 


r 2 u r 

f f 1 (z  - z,'?cos(m;)co.;  [(1+1  , e'3klR 

•'ll  j •'-I.  L R0 


0 


-jk,R, 

+ (1  + jk2RQ)  e L i + 


'•] 


2 (£  - £ * ) cosy 1 


[(£  - £’)2  + pV] 


(3/2) 


2 ir 

f.  L 


2(£-£l)cosy* 


4j  J-it  [(£  - £’)2  + pV] 


(3/2)  d^d^', 


d£d£' 

(A. 21a) 


(A. 21b) 


The  integrand  of  1^  is  non-singular  and  may  be  integrated 
numerically.  Analytical  evaluation  of  1^  yields 


. A£0£ll^tn 

6 P 


{■ 


prr  4-  ^|  (£  - £ j ) 2 + p2if 2' 


* 


pir  + *y<£  - £2)2  + p2ir2 


+ £n 


£ - £. 


£ - £ 


(A. 22) 


which  is  infinite  if  £ = £^  or  £*£2*  These  values  mt.y  oc- 
cur in  the  evaluation  of  (2.25g),  so  has  a non-integrable 
singularity.  The  consequences  of  this  result  are  discussed 
in  Section  A. 9.  We  similarly  express 
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f 


(A. 23) 


U2  - l7  + ^6 

where 


*7  " 


2 TT  /“ 

f f ht-.')p'co.W)c..;  [(l+jklR0)  e‘JklR° 

Am  L R0  ° 


+ (1 + jk2RQ)  e 


_jk2Rol  2p(&  - tMcosy'  i 

J r 2 2 2l(3/2)| 

[tt-4')z  + PV  J 


I 


dCdJl’  . 
(A.  24) 


Thus  U2  also  has  a non-integrable  singularity  and  will  be 
treated  in  Section  A. 9. 

A. 5 Singularity  Analysis  of  the  Integral  Function  u(? 


By  direct  comparison  with  Uq  one  can  immediately  write 


U3  ' J8  + PJ4 


(A. 25) 


where 


2 71  f 

X f [(1+3kiV  e'JklR0 

A*  L Ro 


+ Cl  + jk2RQ)  e * + 


-jk2R0‘ 


r 2 2 2"l^3^2^ 

L(»-r)"  + PVJ 


d£d*'  , 

(A.  26) 
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and  where  1^  is  given  by  (A.19)  . The  integral  Ig  has  a 
non-singular  integrand  and  can  be  integrated  numerically. 

A. 6 Singularity  Analysis  of  the  Integral  Function 

By  analogy  with  one  can  compute 

u4  * + I10  > (A. 27) 

where 


P 1 (p  - P ' )cos (mg) 
>: 


0 


(1  + jkjR0) 


-JkiR0 


+ (1  + jk2RQ)  e 


“^2*0 


] 


2p  (SL  - Z’  )sinY  * 

T 2 2 2*1  ^ ^ ^ 

|U  - A'K  + p £ J 


> d£d£' 

(A. 28) 


I10  - 4siny' 

The  integral 
infinite  if  Ji, 
Section  A. 9. 


£n 


PIT  + y(Z  - ij)2  + p2it2'' 


P7T  + fa  - 1 2)  2 + p2*2' 


+ in 


l-l. 


I - l. 


(a. 29) 


Ig  can  be  computed  numerically , but  I ^ ^ is 
* or  2»  and  is  therefore  considered  in 
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t 


Singularity  Analysis  of  the  Integral  Function  U 


and 


The  integral  function  is  defined  by  (2.26h)  as 


(l 


cos  {mt, ) sin 

„3 


(1  + J^Rq)  e 


-JkiRo 


"j  ^2R0  1 
+ (1  + jk.Rn)  e J 


d?dr 


As  R_  -*•  0 the  integrand  of  Uc  approaches 


2 |U  - V)  + p 5 J 

The  integrand  is  therefore  singular  and  we  compute  as 


U5  " *11  + 1 1 2 


(A. 31) 


where 
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(A. 32a) 
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(3/2)  d^d2' 


--3  |<*  - ^Hn  £p 
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- (A  - &1>&n<£.  - Aj)  - a2-£HnU2 


- 5l)|  • 


(A. 32b) 


The  integral  I ^ is  integrated  numerically  since  the  inte- 
grand is  now  non-singular.  Similarly, 


U5  * *13  + >ll2 


(A. 33) 


where 
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[(1  + ikiV 


-Jk!R0 
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-jk2R0 
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(A. 34) 
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A. 8 Singularity  Analysis  of  the  Integral  Function 

Analysis  of  the  integral  function  is  directly  anal- 
ogous to  that  of  U*J.  We  have 


U6  35  l14%110  ’ 


(A. 35) 


where 
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+ (l  + jk2R0)  e 


^jk2R 


2(2,  - 2,  * ) sinv  * 


[u  - v)1  + pV] 


(3/2) 


f d£d2. ' , 

(A. 36) 


and  where  I^q  is  given  by  (A. 29).  The  integral  I ^ can  be 
evaluated  numerically,  whereas  the  analytical  result  for 
may  be  infinite  and  is  therefore  treated  in  the  next 
section. 

A . 9 Elimination  of  Non-Integrable  Singularities 

The  non-integrable  singularities  which  have  appeared 

P P 

in  the  analyses  of  the  integral  functions  U2,  U2,  U£,  and 
U,  are,  of  course,  also  non-physical.  We  therefore  proceed 


to  determine  how  these  apparent  singularities  may  be  re- 
moved. Consider  now  the  quantity  K which  is  the  sum  of  the 
third  and  fifth  terms  of  (2.25g).  Then  for  a self  term  we 
have 


4ttK 

At 
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-U,cosy  + U-sinv 
6 1 n 2 ' n 

cosy 

-I.-cosy  - I,. — - + Icsiny  + I, siny 
14  ' n lOp  5 ' n 6 ' n 

(A. 37) 


Then  note  that 

cosy 

-I,., — - + I siny 
1 0 p 6 n 


where 


-4cosy  siny 
n 1 n 


0 


4siny  cosy 
n n 


A , 

(A. 38a) 


£n 


pn  +^(£  - S.^2  + p2ir2 


{< 


pi!  + - z2)2  + p2*2 


+ In 


i - a. 


l - 1 


(A. 38b) 


Thus,  the  result  for  K does  converge  as  expected.  The 
numerical  procedure  for  the  evaluation  of  the  integral 
functions  and  is  therefore  to  compute 

U2  - I5  (A. 39a) 
U6  - Iu  . (A. 38b) 
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thus  removing  the  canceling  non-integrable  terms  from  the 
evaluation  of  the  remaining  terms  in  the  integral. 

If  we  next  consider  the  quantity  K which  is  the  sum  of 
the  first  and  third  terms  of  (2.25f),  then  for  a self  term 
we  have,  for  example, 

2K  - Xc(4VYn)<  - X5<Atn,Yn)»5 

Xc(Atn’YnJ  (I9  + I10)  ~ Xs ' Atn,Yn* (I7  + pI6)  » 

(A. 40) 

and 

Xc^Who  * x»<Atn’Yn>f,I6 

■'  (2[At„+1co,Yn+1  + AtncosY11l8lnYm 

- 2[Atn+l8l”Yn+l  + Atn8inYn]co8Y,»)A  • 

(A. 41) 

where  A is  defined  by  (A. 38b)  and  where,  for  a self  term, 
either  m - n or  m«n+l.  By  way  of  example  we  choose  m-n, 
which  gives 

VlO  - Vr6  ' 2Aitn-H!oosYn+l8lnYn  ' 8lnY»+lc°8Y„}  ' 

(A. 42) 

The  trigonometric  functions  in  (A. 42)  do  not  cancel  unless 
Yn  * (which  implies  that  there  is  no  bend  in  the  surface 
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at  t = t ) and  the  result  may  therefore  be  infinite.  The 
n 

difficulty  lies  in  the  mathematical  representation  of  the 
problem  — the  field  point  is  placed  between  two  current 
subdomains  each  of  which  has  an  independent  associated  cur- 


rent  c 

oef  f icient 

and  this  point 

may  be 

a t 

a bend  in 

the 

surf  ac 

e where  Y 

n 

'Vi'  u 16 

relativ 

el  y 

easy 

to  d 
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true 
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in  the 
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would 
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in  the  Cau 

chy  Principal  V 

alue  sense. 

To 

eliminate 

this  p 

roblem,  we 

define  in  this 

situat 

ion 

U5  ■ *7 

(A. 43a) 

< - 1, 

> 

(A. 43b) 

thus  r 

emoving  th 

e singular  term 

s from 

the 

eval 

uatio 

n of  the 

integr 

als.  Whil 

e this  procedur 

e is  pe 

rhap 

s no 

t ent 

irely 

rigorous,  it  has 

been  found  to 

work  sa 

t isf 

acto 

rily. 

* 
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Appendix  B 


DESCRIPTION  AND  LISTING  OF  THE  CODE 

The  computer  code  DBR  is  designed  to  calculate  the  in- 
duced currents  on  a general  body  of  revolution,  which  may  be 
either  a perfectly  conducting  or  a dielectric  body.  Perfectly 
conducting  bodies  may  be  either  '’closed"  or  "open,"  i.e.  the 
body  surface  may  or  may  not  intersectthe  axis  of  the  body  of 
revolution,  respectively.  Dielectric  bodies  must  be  closed. 
The  code  is  not  capable  of  treating  bodies  in  which  the  gen- 
erating arc  forms  a closed  loop  (such  as  for  toroidal  bodies) 
nor  is  it  capable  of  treating  bodies  having  a finite  conduc- 
tivity. The  modifications  necessary  to  include  such  cases, 
however,  are  fairly  straightforward.  The  excitation  which 
induces  current  on  the  body  is  assumed  to  be  a plane  wave 
for  perfect  conductors  or  dielectrics  and/or  a single  delta- 
gap  voltage  source  for  perfect  conductors.  Any  delta-gap 
voltage  source  is  assumed  to  be  (^-independent . 

The  output  of  the  code  for  each  excitation  consists  of 
a description  of  the  input  parameters  and  pertinent  calculated 
parameters,  listings  of  each  Fourier  component  of  current, 
and  listings  of  sums  over  the  computed  Fourier  current  com- 
ponents observed  in  specified  planes  of  constant  <f>. 
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The  code  itself  contains  sufficient  information  in  comment 
cards  for  user  operation.  The  general  purpose  auxiliary 
routines  needed  by  the  code  are  also  listed  in  this  appendix 
for  the  reader’s  convenience. 

B . 1 Program  Operation 

The  program  MAIN  reads  all  input  data  except  for  input 

data  which  describes  the  generating  arc.  The  input  data  is 

adequately  described  by  comment  cards  in  the  program.  The 

subroutine  GENCUR  is  called  to  read  the  data  describing  the 

generating  arc  in  terms  of  the  p-  and  z-coordinates  of  the 

points  t (see  Fig.  2.2).  GENCUR  also  calculates  and  stores 
n 

vectors  corresponding  to  the  quantities  t , , At  , and  Y . 

n+*s  n n 

The  number  of  points  used  should  be  sufficient  to  adequately 
approximate  the  generating  arc  as  well  as  to  represent  the 
variation  of  the  current  on  the  body.  The  body  surface 
should  not  intersect  the  ...vis  of  the  body  of  revolution  ex- 
cept at  the  end  points  of  the  generating  arc.  If  the  surface 
does  intersect  the  axis  at  either  or  both  ends  of  the  gen- 
erating arc,  then  the  points  t^  and/or  tN+^  should  be 
placed  on  the  axis.  The  general,  form  of  the  input  data  for 
a single  case  is  given  on  the  next  page  in  terms  of  the 
program  variables: 
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NFLDS 

LAMBDA 

MODEB 

EPSR 

MODEE 

MUR 

NGQ 

1DB 

THETA(l) 

• 

• 

• 

PHI (1) 

ETHETA ( 1 ) 

EPH 1(1) 

ANTFD ( 1 ) 

THETA(NFLDS)  PHI (NFLDS ) ETHETA (NFLDS ) EPHI(NFLDS)  ANTFD (NFLDS) 
RHO(l)  Z(l) 

RHO ( 2 ) Z ( 2) 

4 

RHO(NPTS)  Z(NPTS) 

9999.0  9999.0 

Note  that  the  number  of  excitations  entered  must  correspond 
to  the  input  parameter  NFLDS.  On  the  other  hand,  the  points 
read  which  describe  the  generating  arc  a^e  terminated  by  the 
presence  of  the  numbers  9999.0  and  the  quantity  NPTS  is  cal- 
culated internally.  Multiple  cases  may  be  run  by  placing  a 
new  data  set  immediately  following  the  first  data  set. 

Once  all  data  is  read  for  a particular  case  MAIN  computes 
the  dimensions  required  to  run  the  case.  If  the  absolute 
dimensions  provided  in  MAIN  are  inadequate , execution  is 
aborted  for  the  case  and  data  is  read  for  a new  case,  if 
present;  otherwise,  SLTN  is  called  to  compute  and  print  the 
results.  The  impedance  matrix  ..nd  drive  vector  are  computed 
for  each  Fourier  component  by  calls  to  the  subroutines  ZMATRX 
and  CVFILL,  respectively.  The  elements  computed  by  these 
routines  are  twice  those  indicated  by  Eq.  (2.25)  and  (2.30). 

The  impedance  matrix  is  inverted  by  the  routine  CSMINV  and 
the  current  coef f io i ent s indicated  in  (2.17)  are  computed  as 


r, 


I = Z-1V  for  each  Fourier  component  by  the  multiplication 
routine  ICRMUL.  The  coefficients  of  the  exponential  series 
(2.17)  are  transformed  to  coefficients  of  a trigonometric 
series  by  the  subroutine  SUMOUT.  The  current  coefficients 
are  then  printed  and  the  process  is  repeated  for  the  next 

l. 

Fourier  component. 

As  mentioned  previously,  the  elements  of  the  impedance 

matrix  calculated  by  the  program  are  twice  that  Indicated 
by  (2.25).  Except  for  this,  the  elements  are  calculated 

using  the  expressions  (2.25)  and  the  duality  relations  (2.14), 
where  appropriate.  The  integral  functions  \p  and  U appearing 
in  (2.25)  are  calculated  by  numerical  integration.  Gaussian 
quadrature  integration  (subroutines  CGQl,  CGQ1T)  is  used  in 
the  t-direction  and  the  order  of  the  approximation  is  controlled 
by  the  input  parameter  NGQ.  Second  order  Gaussian  quad- 
rature (NGQ  = 2)  has  generally  been  found  sufficient.  Gaussian 
quadrature  integration,  however,  represents  a fixed  order 
approximation  and  it  is  not  entirely  satisfactory  for  per- 
forming integrations  in  the  (^-direction  since  the  integrands 
of  the  ’p  and  U functions  may  vary  rapidly  as  a function  of  <J> . 

We  have  therefore  chosen  to  per  form  integrations  in  the  in- 
direction using  adaptive  trapezoidal  rule  integration  (sub- 
routine TRPADP) . The  adaptive  numerical  integration  is 
terminated  when  the  relative  error  between  the  two  most  recent 
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results  is  less  than  1%.  The  adaptive  integration  procedure 
assures  that  the  integrations  in  the  <j>-direction  will  be  per- 
formed accurately  regardless  of  the  body  size  or  Fourier  com- 
ponent considered.  For  self-terms,  in  which  the  integrands  of  the 
ij;  and  U functions  may  be  singular,  the  procedures  indicated 
in  Appendix  A are  followed.  Any  analytically  evaluated  portion 
of  a»  integral  function  U is  included  by  the  function  subrou- 
tine CANALY.  Analytically  evaluated  portions  of  the  ip 
functions  are  included  automatically  in  the  subroutines 

ELLPTC  and  ELLPTR.  For  non-self  terms  the  \p  functions  are 
evaluated  as  indicated  by  Eq.  (A. 13)  in  order  to  increase  the 
convergence  rate  of  the  adaptive  integration  in  the  <{>-direc- 
tion  by  smoothing  the  integrand. 

B . 2 Sample  Case 

A sample  case  is  provided  in  this  section  to  provide  a 
convenient  check  for  user  implementation  only.  The  sample 
run  is  not  intended  to  adequately  model  any  structure,  but 
simply  to  exercise  most  sections  of  the  code.  The  input  data 
for  the  sample  case  follows: 


1 

0 

1.0 

4.0 

90. 

0. 

0.0 

0.0 

.25 

0.0 

.5 

.25 

.5 

.5 

.3 

.65 

.15 

.65 

0.0 

. 6 

9999. 

9999 

1.0 

-1.0  , 0.0  0.0 


The  resulting  output  of  the  program  is  presented  on  the  fol- 
lowing pages. 


IBIS  PACE  IS  BfcST  QUALITY  PRACTICABLI 
I3WA  OWY  JUHKSHJSD  TO  DDC  


*************************** 


* * 

» CASE  NUMBER:  1 * 

* * 

* EXCITATION  NUMBER:  1 * 

* * 


POINT 

RHO 

2 

DELTA-T 

GAMMA 

1 

0 .0  00000E+00 

0.000000E+00 

1.5 70796 E+00 

1.250000E-01 

0 , 0 00000E+00 

2.500000E-01 

2 

2.500000E-01 

0,0  00000E+  00 

7.853982E-01 

3.750000E-01 

1.250000E-01 

3.535534E-01 

3 

5.000000E-01 

2. 5000 00 E- 01 

5.000000E-01 

J.750000E-01 

2.500000E-01 

0.0  00000E+00 

4 

5.000000E-01 

5.000000E-01 

-9.272952E-01 

4,000000 E~ 01 

5. 7 500 00 E- 01 

2.500020E-01 

5 

3.000000E-01 

6. 5000  00  E-*  01 

-1.570796E+00 

2.250000E-01 

6.500000E-01 

1 .5  00000E-01 

6 

1.500000E-01 

6 . 5 00000E-01 

-1.892547E+00 

7 .500000E-02 

6.2 J0000E-01 

1.581139E-01 

7 

0.000000E+00 

6.000000E-01 
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5Kf  BBST  QUAim  PiucriciiM 

**»**..*.**«**  WJ*  OQKf  JIBQUSHKD  TO  HOC 

* * 

* INPUT  DATA  * 

* * 

************** 

MAXIMUM  NUMBER  OF  MODES  TO  BE  USED  « 2 

ORDER  OF  GAUSSIAN  QUADRATURE  - 2 

FREE  SPACE  WAVELENGTH  - 1.0000000E+00 

A DIELECTRIC  BODY  HAS  BEEN  ASSUMED: 

RELATIVE  DIELECTRIC  CONSTANT  " 4.0000000E+00 

RELATIVE  PERMEABILITY  - 1.0000000E+00 

INCIDENT  FIELD  DATA: 

THETA  * 9 .0 000000E+01 

PHI  - 0 .0  000000E+00 

E-THETA  - -1. 0000000 E+ 00 

E-PHI  - 0.0000000E+00 

* * 

* COMPUTED  DATA  * 

* * 

***************** 

PARAMETERS  OF  FREE  SPACE: 

DIELECTRIC  CONSTANT  - 8 .8541850E-12 

PERMEABILITY  - 1 .2 56637 0E- 06 

WAVENUMBER  • 6 .28318 50E+00 

SPEED  OF  LIGHT  ■ 2 .9 979250E+ 08 

PARAMETERS  OF  THE  DIELECTRIC  BODY: 

DIELECTRIC  CONSTANT  - 3 .5 416740E- 11 

PERMEABILITY  - 1 . 2 566370 E- 06 

WAVENUMBER  - 1 .2 566370 E+01 

SPEED  OF  LIGHT  - 1 .4 989630E+08 

WAVELENGTH  « 5 .0 000000E- 01 

OMEGA  ~ 1.8836520E+09 

FREQUENCY  - 2 .997925 0C* 08 
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"T"  COMPONENT  OF  ELECTRIC  CURRENT 
FOR  THETA  DIRECTED  PORTION  OF  INCIDENT  FIELD 
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B.3  Program  Listing 


c 

PROGRAM  DRR 

BR- 

10 

c 

BR- 

20 

c 

VERSION  4.01 

BR- 

30 

c 

BK- 

4V 

0 

JULY  20,  1978 

BR- 

50 

(j***  * ***************  ************  ******  ****  *********  ****  ****  *******  **  **  **RR— 

60 

c 

BR- 

70 

c 

PROGRAM  DBR  CALCULATES  THE  CURRENTS  (TOTAL  AND  DENSITY)  FOR  A 

BR- 

80 

c 

BODY  OF  REVOLUTION,  EITHER  DIELECTRIC  OR  PERFECTLY  CONDUCTING. 

BR- 

90 

c 

BR- 

100 

c 

BR- 

110 

c 

PROGRAM  "MAIN.'  READS  ALL  THE  INPUT  DATA  EXCEPT  FOR  DATA 

BR- 

120 

c 

CONCERNING  THE  GENERATING  CURVE  OF  THE  BODY  OF  REVOLUTION. 

BR- 

130 

c 

SUBROUTINE  "GENCUR"  IS  CALLED  TO  READ  THIS  DATA. 

"MAIN."  THEN 

BR- 

140 

c 

COMPUTES  THE  DIMENSIONS  NECESSARY  TO  RUN  THE  PROGRAM  FOR  THE  CASEBR- 

150 

c 

IN  QUESTION  AND  ABORTS  EXECUTION  OF  THE  CASE  IF  THE  ABSOLUTE 

BR- 

160 

c 

DIMENSIONS  PROVIDED  IN  "MAIN."  ARE  INSUFFICIENT. 

IF  THE 

BR- 

170 

c 

DIMENSIONS  ARE  ADEQUATE  "MAIN.*  CALLS  SUBROUTINE 

"SLTN"  TO 

BR- 

180 

c 

SOLVE  THE  PROBLEM. 

BR- 

190 

c 

BR- 

200 

c 

AUXILIARY  ROUTINES  NEEDED: 

BR- 

210 

c 

BR- 

220 

c 

CGQ1  - INTEGRATION  OF  A COMPLEX  FUNCTION  BY 

GAUSSIAN  QUAD. 

BR- 

230 

c 

TRPADP  - INTEGRATION  OF  A COMPLEX  FUNCTION  BY 

ADAPTIVE 

BR- 

240 

c 

TRAPEZOIDAL  RULE 

BR- 

250 

c 

CSMINV  - INVERSION  OF  A COMPLEX  MATRIX 

BR- 

260 

c 

ICRMUL  - COMPLEX  MATRIX  MULTIPLICATION 

BR- 

2 70 

0 

BESSEL  - COMPUTATION  OF  BESSEL  FUNCTION  OF  ARBITRARY  ORDER 

BR- 

280 

c 

ELIC1K  - COMPUTES  COMPLETE  ELLIPTIC  INTEGRAL 

OF  FIRST  KIND 

BR~ 

290 

c 

BR- 

3 00 

c 

ASSIGNMENTS: 

BR- 

310 

c 

BR- 

3 20 

c 

READ  FROM  DSK:  *2 

BR- 

330 

c 

WRITE  PRINTED  OUTPUT  ON  LPT:  *5  THROUGH  *15 

BR- 

340 

c 

(AS  NECESSARY:  FOR  MULTIPLE  EXCITATIONS) 

BR- 

350 

c 

BR- 

360 

c 

BR- 

370 

£********************»«********  *********************************  ******  * * BR“* 

380 

IMPLICIT  COMPLEX  (C) 

BR- 

39  0 

REAL  LAMBDA, MU, MUR 

BR- 

400 

c 

-LOCATE  1" 

BR- 

410 

DIMENSION  C(10296) ,R(306) 

BR- 

4 20 

COMMON/M  BE/MODEB , MODE  E 

BR- 

430 

COMMON/PNT/NPTS ,NUNKT,NUNKPH,NUNK2 

BR- 

440 

CQMMON/INT/NMODEM, IDB ,NGQ 

BR- 

450 

COMMON/WVL/LAMBDA,EPSR,MUR 

BR- 

460 

COMMON/CAS/ICASE 

BR- 

470 

c 

"LOCATE  2" 

BR- 

4 80 

NDMSNC=10296 

BR- 

490 

NDMSNR“306 

BR- 

508 

c 

BR- 

510 
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-NOTE:  WHEN  THE  DIMENSION  OF  "C"  OR  "R"  IS  CHANGED,  THE  VALUE  OF 
"NDMSNC"  OR  "NDMSNR" , RESPECTIVELY , MUST  BE  CHANGED  TO  THE 
APPROPRIATE  VALUE  FOR  PROGRAM  OPERATION. 

INPUT  DATA  DESCRIPTION: 

{ALL  DATA  IS  READ  IN  "FREE"  FORMAT) 

NFLDS  = NUMBER  OF  DIFFERENT  EXCITATIONS  TO  BE  SOLVED  FOR 
THIS  STRUCTURE. 

MODEB  = STARTING  FOURIER  COMPONENT  NUMBER 
MODEE  “ ENDING  FOURIER  COMPONENT  NUMBER 

IDB  = A CONTROL  VARIABLE: 

* 1,  IF  THE  STRUCTURE  IS  A DIELECTRIC  BODY 

* 0,  IF  THS  STRUCTURE  IS  A PERFECT  CONDUCTOR 

ANTFD  “ NUMBER  OF  SUBDOMAIN  WHICH  IS  FED  BY  A DELTA  GAP  SOURCE. 
IF  =0,  THERE  IS  NO  SOURCE.  THE  SECOND  COORDINATE 
POINT  READ  IN  IS  CONSIDERED  TO  BE  THE  FIRST  SUBDOMAIN, 
AND  SO  ON. 

EPSR  * RELATIVE  DIELECTRIC  CONSTANT  OF  THE  BODY 
(IF  DIELECTRIC) 

MUR  = RELATIVE  PERMEABILITY  OF  THE  BODY  (IF  DIELECTRIC) 
LAMBDA  “ WAVELENGTH  (IN  METERS)  IN  FREE  SPACE  . 

NGQ  * ORDER  OF  GAUSSIAN  QUADRATURE  INTEGRATION  TO  BE  USED 
THETA  “ SPHERICAL  COORDINATE  ANGLE  OF  INCIDENCE  (IN  DEGREES) 

PHI  = SPHERICAL  COORDINATE  ANGLE  OF  INCIDENCE  (IN  DEGREES) 
ETHETA  “ SIGNED  MAGNITUDE  OF  THE  E-THETA  INCIDENT  FIELD 
EPHI  * SIGNED  MAGNITUDE  OF  THE  E-PHI  INCIDENT  FIELD 


DO  40  ICASE  - 1,  50 

READ(2 , 10000, END=  50) NFLDS , MODEB , MODE  E, NGQ, IDB 
READ(2, 10001)  LAMBDA,  EPSR,  MUR 
NMODEM=MODEE-MODEB+l 
. I2-NFLDS 
I3=NFLDS+I2 
I4*NFLDS+I3 
I5“NFLDS+I 4 
I6“NFLDS+I5 
DO  10  I * 1,  NFLDS 

READ(2, 10001)  THETA,  PHI,  ETHETA,  EPHI,  ANTFD 
R ( I) “THETA 
R ( 1 2+1 ) “PHI 
R(I3+I) “ETHETA 
R(I4+I) “EPHI 
RU5+I)  “ANTFD 
10  CONTINUE 

CALL  GENCUR l R ( 16  + 1) , 1 6 , NDMSNR , N PT  S , N REQR ) 
NREQR=N  REQR+1 6 
NUNKT=N  PTS- 2 


8R- 

520 

BR- 

530 

BR- 

540 

BR- 

550 

BR- 

560 

BR- 

570 

BR- 

580 

BR- 

590 

BR- 

600 

BR- 

610 

BR- 

620 

BR- 

630 

BR- 

640 

BR- 

650 

BR- 

660 

BR- 

670 

BR- 

680 

BR- 

690 

BR- 

70  0 

BR- 

710 

BR- 

720 

BR- 

730 

BR- 

74  0 

BR- 

750 

BR- 

/60 

BR- 

770 

BR- 

780 

BP- 

790 

BR- 

800 

BR- 

810 

BR- 

820 

BR- 

83  0 

BR- 

840 

BR- 

850 

BR- 

860 

BR- 

870 

BR- 

880 

BR- 

890 

BR- 

90  0 

BR- 

910 

BR- 

920 

BR- 

930 

BR- 

940 

BR- 

950 

BR- 

960 

BR- 

970 

BR- 

980 

BR- 

990 

BR- 

1000 

BR- 

1010 

BR- 

1020 

BP.- 

1030 

BR- 

1040 

BR- 

1050 

BR- 

1060 

BR- 

1070 

BR- 

1080 

tsR~ 

1090 

BR-- 

1100 

BR- 

1110 

BR- 

1120 
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NUNKPH=NPTS-1  BR- 

IF(IDB  .EQ.  0)  GO  TO  20  BR- 

NUNKT=>NUNKT*2  BR- 

NUNKPH=NUNKPK*2  BK- 

20  NUNK“NUNKT+NUNK  PH  BR- 

NUNK2=NUNK*NUNK  BR- 

NREQC=NUNK*NFLDS*6+NUNK2  BR- 

' IF(  NREQC  .GT.  NDMSNC  .OR.  NREQR  .GT.  NDMSNR)  GO  TO  30  BR- 

NPTSM1=NPTS-1  BR- 

MR2*NFLDS+1  BR- 

M R3  *N  FLDS+MR2  BR- 

MR4  *N  FLDS+MR3  BR- 

MR5*NFLDS+MR4  BR- 

MR6-NFLDS+MR5  BR- 

MR7 “NPTSM1+MR6  BR- 

MR8-NPTSM1+MR7  BR- 

M R9  *N  PTS  Ml  +M  R8  BR- 

MR10-NPTSM1+MR9  BR- 

MR11“NPTS+MR10  BR- 

C MR12»NPTS+MRll  BR- 

MC2-NUNK2+1  BR- 

MC3«NUNK*NFLDS+MC2  BR- 

MC4«NUNK*NFLDS+MC3  BR- 

MC5-NUNK*NFLDS+MC4  BR- 

MC6«NUNK*NFLDS+MC5  BR- 

C MC7«NUNK»NFLDS+MC6  BR-* 

CALL  ZERO(C, NREQC)  BR- 

CALL  SLTN{  3R- 

5 C(l)  ,C(MC2) ,C(MC3),C(MC4),C(MC5),C(MC6),  BR- 

5 R (1 ) , R (MR2) , R(MR3 ) , R(MR4) , R(MR5 ) , R(MR6 ) , R(MR7 ) , R(MR8 ) , BR- 

$ R(MR9 ) , R(MR10) , R(MRll) ,NUNK ,NFLDS  ) BR- 

GO  TO  40  BR- 

30  CONTINUE  Bk- 

WKlTtl(5, 10002)  ICASE,  NDMSNC,  NDMSNR,  NREQC,  NREQR  BR- 

40  CONTINUE  BR- 

50  CONTINUE  BK- 

10000  FORMAT (91)  BR- 

10001  FORMAT (9E)  BR- 

10002  FORMAT ( ' CASE  NUMBER' 13,'  EXECUTION  ABORTED:  DIMENSIONS  INSUFFICIBR- 
IENT' // ' DIMENSION  GIVEN  IN  PROGRAM: ’ //5X, ' C ( '16 , ’ ) *5X, ’ R( '1 6, ' ) '//BR- 
#'  DIMENSIONS  REQUIRED  FOR  THIS  CASE : ' //5X, 'C ( '16 , ' ) '5X, ' R( '1 6, ' ) ’ )BR- 

STOP  BR- 

END  8k- 

SUBROUTINE  GENCUR( S ,NFLDS4 , NDMSNS ,NPTS, NREQR)  BR- 

»*« »r***»********«********************************«*««**************** 

BR- 


SUBROUTINE  "GENCUR"  READS  THE  GENERATING  CURVE  DATA  FOR  THE  BODY  BR- 
OF  REVOLUTION  AND  USES  THIS  INFORMATION  TO  COMPUTE  THE  DIMENSIONS  BR- 
NECESSARY  OF  THE  REAL  VECTOR  "S".  IT  ALSO  USES  THE  GENERATING  BR- 
CURVE  DATA  TO  COMPUTE  INTERMEDIATE  COORDINATE  LOCATIONS,  VALUES  BR- 
OF  STEP  SIZE,  AND  ANGLES  CORRESPONDING  TO  THE  DIRECTION  OF  EACH  BR- 


ELEhENTAL  SURFACE.  BR- 

BR- 

*********  ft********************  * * **  *********  **  ****  ****  **  ***  **  ****  **  ****  * g 

DIMENSION  S(l)  BR- 

NDMSN=NDMSNS-NFLDS4  BR- 

BR- 

INPUT  DATA  DESCRIPTION:  BR- 

BR- 

RHO  » CYLINDRICAL  COORDINATE  RHO  (IN  METERS)  BR- 

BR- 


1130 
1140 
1150 
1160 
1170 
1180 
1190 
1200 
1210 
1220 
1230 
1240 
1250 
1260 
1270 
12U0 
129  0 
1300 
1310 
1320 
1330 
1340 
1350 
1360 
1370 
1380 
1390 
1400 
1410 
1420 
1430 
1440 
1430 
1460 
1470 
1480 
1490 
1500 
1510 
1520 
1530 
1540 
1330 
1560 
1570 
1580 
1590 
1600 
1610 
1620 
1630 
1640 
1650 
1660 
1670 
1680 
1690 
1700 
1710 
1720 
1730 
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Z = CYLINDRICAL  COORDINATE  Z (IN  METERS) 

NOTE:  THE  LAST  RECORD  OP  THE  GENERATING  CURVE  DATA  MUST  BE 
9999.  9999. 

(THE  SPACING  IS  ARBITRARY). 


1*1 

DO  10  K * 1,  NDMSN 
10  S (K) *0,0 

20  READ(2 (1 0000)  RHO,  Z 
NPTS“I-1 

IF(RHO  ,EQ.  9999.  .AND.  Z .EQ.  9999.)  GO  TO  30 
K*2*I 

IF(K  .GT.  NDMSN)  GO  TO  20 
S(K-l) *RHO 
S(K)*Z 
1*1+1 
GO  TO  20 

COMPUTE  THE  DIMENSION  OF  "S*  REQUIRED  FOR  THIS  CASE 

30  NREQR«6*NPTS-4 

IF (NREQR  .GT.  NDMSN)  RETURN 
NSTRT*  4*NPTS- 4 

REORDER  THE  INPUT  DATA 

DO  40  1*1,  NPTS 
K-2*I 

NI*NSTRT+I 
S (NI) *S (K-l) 

S (NPTS+NI) *S (K) 

40  CONTT‘JUE 

NPTSM1-NPTS-1 
DO  50  1*1,  NPTSMl 
NI-NSTRT+I 
RI*S (NI) 

RIPl*S (NI+1) 

NI*NPTS+N I 
2 I*S (NI) 

Z IPl*S (NI+1) 

STORE  RHO  AND  Z HALFWAY  POINTS 

S (I)  * (RI+RIPl) /2.0 
S ( I+NPTSM1 ) * (Z I+ZIP1) /2 ,0 
DIIPi*SQRT( (RIP1-RI) **2+(ZIPl~ZI) **2) 

ARGN* (RIPl-RI) /DIIP1 
ARGD* (ZIPl-ZI) /DIIP1 

STORE  DELTA-T  VECTOR 

S (I+2wNPTSMi ) *D IIP1 

STORE  GAMMA  VECTOR 

S(I+3*NPTSM1) =ATAN2(ARGN,ARGD) 

50  CONTINUE 
10000  FORMAT (9E) 

RETURN 


BR-  1740 
BR-  1750 
BR-  1760 
BR-  1770 
BR-  1780 
BR-  17  90 
BR-  1800 
BR-  1810 
BR-  1820 
BR-  1830 
BR-  1840 
BR-  1850 
BR-  1860 
BR-  1870 
BR-  1880 
BR-  1890 
BR-  1900 
BR-  1910 
BR-  1920 
BR-  1930 
BR-  19  40 
BR-  1950 
BR-  19  60 
BR-  1970 
BR-  19  80 
BR-  1990 
BR-  2000 
BR-  2010 
BR-  2020 
BR-  2030 
BR-  2040 
BR-  2050 
BR-  20b0 
BR-  20/0 
BR-  2080 
BR-  2090 
BR-  2100 
BR-  2110 
BR-  2120 
BR-  2130 
BR-  2140 
BR-  2150 
BR-  216  0 
BR-  217  0 
BR-  2180 
BR-  2190 
BR-  2200 
BR-  2210 
BR-  2220 
BR-  2230 
BR-  2240 
BR-  2250 
BK-  2260 
BR-  2270 
BR-  228  0 
BR-  2290 
BR-  2300 
BR-  2310 
BR-  2320 
BR-  2330 
BR-  2340 
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END  BR- 

SU8ROUTINE  SLTN(  BR- 

$ CZ,CIT,CIP,CVT,CVP,CITOLD,  BR- 

$ THETA, PHI , ETHETA,EPHI ,ANTFD, RHOPH , ZPH,DELTAT,  BR- 

$ GAMMA,RHO,Z,NUNK ,NFLDS  ) BR- 

*****************************************  ***************************** *BR— 

BR- 

SUBROUTINE  "SLTN"  COMPUTES  AND  OUTPUTS  THE  RESULTS.  ' BR- 

PRINTED  RESULTS  ARE  OUTPUT  TO  DEVICES  NUMBER  5 THROUGH  15.  BR- 

(FOR  MULTIPLE  EXCITATIONS,  EACH  CASE  IS  OUTPUT  TO  A NEW  BR- 

DEVICE  #,  BEGINNING  WITH  DEVICE  #5)  BR- 

BR- 

************************************************  ********************* **gj^_ 

IMPLICIT  COMPLEX  (C)  BR- 

REAL  LAMBDA, MUR, MU1, MU 2 BR- 

DIMENSION  POL(2) ,FTC(2)  BR- 

DIMENSION  CZ(NUNK,NUNK) ,CIT(NUNK,NFLDS ) , CIP (NUNK,NFLD6)  BR- 

DIMENSION  CVT(NUNK,NFLDS) ,CVP{NUNK,NFLDS)  BR- 

DIMENSION  CITOLD(NUNK,NFLDS, 2)  BR- 

DIMENSION  THETA (1 ) ,PHI (1 ) , ETHETA (1 ) , EPHI (1 ) , ANTFD (1)  ' BR- 

DIMENSION  RHO (1 ) , Z (1 ) , RHOPH (1 ) , ZPH (1 ) , DELTAT(1 ) ,GAMMA(1 ) BR- 

COMMON/MBE/MODEB,MODEE  BR» 

COMMON/SCAFAC/ISCALE  BR- 

COMMON/PNT/NPTS ,NUNKT ,NUNKPH,NUNK2  BR- 

COMMON/ INT/NMODEM ,IDB ,NGQ  BR- 

COMMON/WV  L/ LAMBDA , EPS  R, MUR  BR- 

COMMON/CAS/ICASE  BR- 

COMMON/FRQ/P I , AKl , AK  2 , SLl , SL2 , OMEGA  BR- 

COMMON/PRM/EPSl,EPS2,MUl ,MU2  BR- 

DATA  POL/ 'THETA' , ' PHI  '/  BR- 

DATA  FTC/' SIN',' COS',  BR- 

NPTSM1 “NPTS-1  BR- 

ZER=*0,0  BR- 

PI»3. 1415926536  BR- 

AKl « 2. 0 *P I/LAMBDA  BR- 

SL1*2.997925E8  BR- 

0MEGA«AK1*SL1  BR- 

FREQ-OMEGA/2.0/PI  BR- 

KU1-PIM.0E-7  BR- 

EPS1*1.0/ (MU1*SL1*SL1)  BR- 

IF (I DB  .EQ.  0)  GO  TO  10  BR- 

EPS2»EPS1*EPSR  BR- 

MU2»MU1*MUR  BR- 

SL2=1 .0/SQRT{MU2*EPS2)  BR- 

AK2=OMEGA/SL2  BR- 

WVL2«2.0*PI/AK2  BR- 

10  CONTINUE  BR- 

BR- 

PARAMETER  PRINT-OUT  OPTION  BR- 

BR- 

GO  TO  50  BR- 

BR- 

DO  40  JJ*1,NFLDS  BR- 

JW-JJ+4  BK- 

WRITE(JW, 10036)  ICASE,  JJ  BR- 

WRITE(JW, 10000)  BR- 

DO  30  J=1,NPTS  BR- 

WRITE(JW, 10001)  J,  RHO(J) , Z(J)  BR- 

IF( J .EQ.  NPTS)  GO  TO  20  BR- 

WRITE(JW, 10002)  RHOPH (J) , ZPH(J) , DELTAT(J),  GAMMA (J)  BR- 

20  CONTINUE  BP.- 


2350 

2360 

2370 

2380 

2390 

2400 

2410 

2420 

2430 

2440 

2450 

2460 

2470 

2480 

2490 

2500 

2510 

2520 

2530 

2540 

2550 

2560 

2570 

2580 

2590 

2600 

2610 

2620 

2630 

2640 

2650 

2660 

2670 

2680 

2690 

2700 

2710 

2720 

2730 

2740 

2750 

2760 

2770 

2780 

2790 

2800 

2810 

2820 

2830 

2840 

2850 

2860 

2870 

2880 

2890 

2900 

2910 

2920 

2930 

2940 

2950 
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IF(NFLDS  .GT.  11)  WRlTE(JW,l  0033) 

BR-  2960 

30 

CONTINUE 

BK-  29  IV 

40 

CONTINUE 

BR-  2980 

50 

CONTINUE 

BR-  2990 

DO  66  J=*1,NFLDS 

BR-  3000 

JW»J+4 

BR-  3010 

WRITE (JW, 10003) 

BR-  3020 

WRITE (JW, 10004)  NMODEM,  NGQ 

BR-  3030 

WRITE (JW ,10005)  LAMBDA 

BR-  3040 

IF(IDB  .EQ.  0)  WRITE(JW, 10006) 

BR-  3050 

IF(IDB  .EQ.  1)  WRITE (JW ,1 0007)  EPSR,  MUR 

BR-  3060 

WRITE(JW ,10008)  THETA (J) , PHI(J),  ETHETA(J) , EPHI(J) 

BR-  3070 

IANTFD*I FIX( ANTFD (J) ) 

BR-  3080 

IF(IANTFD  .NE.  0)  WRITE (JW ,1 0025)  RHO(IANTFD+l) , 2 (IANTFD+1) 

BR-  3090 

WRITE(JW, 10017) 

BR-  3100 

WRITE(JW, 10022) 

BR-  3110 

WRITE(JW, 10018)  EPSl , MU1,  AKl,  SL1 

BR-  3120 

IF(IDB  .EQ.  1)  WRITE (JW ,10019) 

BR-  3130 

IF(IDB  .EQ.  1)  WR1TE(JW ,10018)  EPS2 , MU2,  AK2,  SL2 

BR-  3140 

IF ( IDB  .EQ.  1)  WRITS(JW ,10023)  WVL2 

BR-  3150 

WF.ITE(JW, 10020)  OMEGA,  FREQ 

BR-  3160 

60 

CONTINUE 

BR-  3170 

M0ONLY=1 

bR-  3180 

MlONLY“l 

BR-  3190 

DO  70  JJ=*1  ,NFLDS 

BR-  3200 

IF(EPHI(JJ)  .NE.  0.0  .OR.  ETHETA(JJ)  .NE.  0.0)  M0ONLY-0 

BR-  3210 

70 

IF(THETA(JJ)  .NE.  0.0  .AND.  THETA (J J)  . NE.  180.0)  M1ONLY-0 

BR-  3220 

IF(M10NLY  .EQ.  1)  NMODEM-) 

BR-  3230 

IF (M0ONLY  .EQ.  1)  NMODEH-1 

BR-  3240 

DO  340  JM-1, NMODFto 

BR-  3250 

MODE=JM-l+MODEB 

BR-  326  0 

IF(M10NLY  .EQ.  1)  MODE-1 

BR-  3270 

I F (M0ONLY  .EQ.  1)  MODE-0 

BR-  3280 

I F ( JM  .NE.  1)  CALL  ZERO < CZ ,NUNK*NUNK) 

BR-  3290 

CALL  ZMATRX (C Z ,RHO , Z , RHOPH , ZPH , DELTAT , GAMMA , NUNK , MODE) 

BR-  3300 

C 

BR-  3310 

C777? ???????? 7? ?????? 7 7 7? 7? 7? ?????? 77? 7? 7? 7?  ?????????????????????? V?????BR-  3320 

C 

MATRIX  PRINTOUT  OPTION 

BR-  3330 

C???????????7???????7???7??7??7? 7777? ?????? 7? 7777? 77 7? 777? 77 ??????? 77 777BH-  3340 

C 

BR-  3350 

C 

WRITE (5, 11200)  ( ( (II,JJ,CZ(II,JJ) ) ,JJ-1, NUNK) ,11-1- NUNK) 

BR-  3360 

c 

BR-  3370 

CALL  CVFILL(CVT, CVP,RHO,Z, RHOPH, ZPH, DELTAT, GAMMA, THETA, PHI, 

BR-  3380 

1 ETHETA, EPHI.ANTFD, NUNK ,NFLDS, MODE) 

BR-  3390 

c 

BR-  3400 

C???????????????????????????????????????????????????????????????????????BR-  3410 

C 

DRIVE  MATRIX  PRINTOUT  OPTION 

BR-  342* 

C??????? ????????????????????????????????? 77? ??????? 7? 7? 7? 7? 7? 7? 7? 77? 7? 77BR-  3430 

C 

BR-  3440 

C 

WRITE(5, 11100)  ( (CVT( I ,1 ) ,CVP(I,1) ), 1-1, NUNK) 

BR-  3450 

C 

BR-  3460 

I SCALE- 0 

BR-  3470 

CALL  CSMINV(CZ, NUNK, NUNK , CDETRM ,ACOND , IERR) 

BR-  3480 

c 

BR-  3490 

C77? 777? 77 7? ??????????????????? 7? 77 7? 77 77 7? 7 77 7 7? 7? 7? 7? 77? 7? 7? 7? 7? 7? 7777BR-  3530 

C 

INVERSE  MATRIX  PRINTOUT  OPTION 

BR-  3510 

C777777777 ???????????? ??????????? ????? 7? ?????????? 77777 777? ??????? 7? 7? 77BR-  3320 

C 

BR-  3530 

C 

WRITE (5, 11 200)  ( ( (II,JJ,CZ(II,JJ) ),JJ-1, NUNK) ,11-1, NUNK) 

BR-  3540 

C 

BR-  3550 

POWER-~FLOAT ( ISCALE) *10. 

BR-  3560 
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ITHETA-1 

CALL  ICRMUL(CZ,CVT/  CIT,NUNK,NUNK  ,NUNK  , NFL  OS) 

CALL  SUMOUT(CIT,NUNK . NFL DS , MODE, I THETA) 

ITHETA*0 

CALL  ICRMUL(CZ,CVr,CIP,NUNK, NUNK, NUNK, NPLDS) 

CAIL  SUMOUT(CIP,nUNK,NFLDS,MODE,ITHETA) 

80  CONTINUE 

CZERO“CMPLX  (0.0 ,0,0) 

NMODEl*N MODEM- 1 

CURRENT  NORMALIZATION  OPTION 

CALL  NORMAL ( C IT , ETHETA , E PH I , NUNK , N PL DS ) 

CALL  NORMAL ( Cl P, ETHETA , EPH I , NUNK  ,NFL DS ) 

90  CONTINUE 
NTE-NUNKT 

IF(IDB  .EQ.  1)  NTE-NTE/2 

NPHE-NTE+1 

DO  330  J=*l , NFLDS 

JW-J+4 

M0  PR" 0 

M1PR-0 

IF(EPHI(J)  .EQ.  0.  .AND.  ETHETA(J)  . EQ.  0.)  M0PR-1 

IF(THETA(J)  .EQ.  0.  .OR.  THETA { J ) .EQ.  180.)  MlPR-1 

IF (M0 PR  .EQ.  1)  M1PR-0 

IF (M0 PR  .EQ.  1)  GO  TO  110 

IF (Ml  PR  .EQ.  1)  GO  TO  110 

DO  100  JSUM=*1,  NUNK 

CVT(JSUM,J)»CITOLD(JSUM ,J,1) 

CVP(JSUM,J)=CITOLD(JSUM,J,2) 

OBI -1.0 

OB2 -S IN i FLOAT (MODE) *3.1415926/2.0) 

IF (MODE  .EQ.  0)  OB2-1.0 
OBSERV-OEl 

IF(JSUM  .GT.  NTE)  OBSERV-OB2 
IF (JSUM  .GT.  2*NTE+NPHE)  OBSERV-OBl 

Cl TOLD ( JSUM  , J ,1 ) -C ITOLD ( J SUM , J , 1 ) +CIT(JSUM,J) *CMPLX (OBSERV,0 . ) 
OBSERV-OB2 

IF(JSUM  .GT.  NTE)  OBSERV-OBl 
IF (JSUM  .GT.  2“NTE+NPHE)  OBSERV-OB2 

CITOLD(JSUM,J,2) *CITOLD(JSUM ,J,2) +CIP(JSUM,J) *CMPLX (OBSERV,0. ) 
100  CONTINUE 

IF(MODE  .EQ.  0)  GO  TO  110 
LOC2-NTE+1 

CALL  RMS ERR (C ITOLD (1 , J ,1 ) ,CVT(1,J) ,NTE,ETTPC) 

CALL  RMSERR(CITOLB(LOC2,J,l) ,CVT(LOC2,J) ,NPHE,EPTPC) 

CALL  HMSERK(CITOLD(l,J,2) ,CVP(1,J) ,NTE,ETPPC) 

CALL  KMSERR(CITOLD(LOC2,J,2) ,CVP(LOC2,J) ,NPHE,EPPPC) 

IF (IDB  .PO.  0)  GO  TO  110 
LOC1-LOC2V  'Hi’ 

LOC2-LOC1+N'.  ' 

CALL  RMSERR(C ITOLD ( LOCI , J , 1 ) , CVTfLOCl , J)  , NTE,  6TTDI) 

CALL  RMSERR(C ITOLD (i.OC2  , J , 1 ) ,CVT(LOC2,J)  ,NPHE,EPTDI) 

CALL  RMSERR (Cl TOLD (LOCI ,J ,2 ) , CVP(LOCl , J) ,NTE, ETPDI) 

CALL  RMSERR(CITOLD(LOC2,J,2) , CVP(LOC2 ,J0 ,NPHE,EPPDI) 

110  CONTINUE 
IFINAL-0 
120  CONTINUE 

IF(M0PR  .EQ.  1 .AND.  MODE  .NE.  0)  GO  TO  320 
IF(M1PR  .EQ.  1 .AND.  MODE  .NE.  1)  GO  TO  320 


BR-  3570 
BR-  3580 
BR-  3590 
BR-  3600 
BR-  3610 
BR-  3620 
BR-  3630 
BR-  3640 
BR-  3650 
Bk-  job» 
BR-  3670 
BR-  3680 
BR-  3690 
BR-  3700 
BR-  3710 
BR-  3720 
BR-  373  0 
BR-  3740 
BR-  3750 
BR-  3760 
BR-  3770 
BR-  3780 
BR-  3790 
BR-  3800 
BR-  3810 
BR-  3820 
BR-  3830 
BR-  3840 
BR-  3850 
BR-  3860 
BR-  3870 
BR-  3880 
BR-  3890 
BR-  3900 
BR-  3910 
BR-  3920 
BR-  3930 
BR-  3940 
BR-  3950 
BR-  3960 
BR-  3970 
BR-  3980 
BR-  3990 
BR-  4000 
BR-  4010 
BR-  4020 
BK-  4030 
BK-  4040 
BR-  4050 
BR-  4060 
BR-  4070 
BR-  4080 
BR-  4090 
BR-  4100 
BR-  4110 
BR-  4120 
BR-  4130 
BR-  4140 
BR-  4150 
BR-  416  0 
BR-  4170 
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DO  310  ITP=1,2  BR- 

IF(M0PR  .EQ.  lj  GO  TO  130  BR- 

IPlITP  .EQ.  1 .AMD.  E THETA (J)  .EQ.  0.)  GO  TO  310  BR- 

130  IF(ITP  .F.Q.  2 .AND.  EPHI(J)  .EQ.  0.)  GO  TO  310  BR- 

PSENS=POL( ITP)  BR- 

IF(ITP  » SQ.  1)  GO  TO  150  BR- 

DO  140  ITR=1,NUNK  BR- 

CIT(ITR,J)=CIP(ITR,J)  BR- 

140  CONTINUE  BR- 

150  CONTINUE  BR- 

PHI1-0.  BR- 

PHI2-90.  BR- 

ET1-ETTPC  BR- 

ET2-ETTDI  BR~ 

EP1-EPTPC  BR- 

EP2-EPTDI  BR- 

1F(ITP  .EQ.  1)  GO  TO  160  BR- 

ET1-ETPPC  BR- 

EP1-EPPPC  BR- 

ET2-ETPDI  ' BR- 

EP2-EPPDI  BR- 

160  CONTINUE  BR“ 

IF (IFINAL  ,NE.  1)  GO  TO  180  BR- 

DO  170  1TR“ 1 ,NUNK  BR- 

170  CIT\ITR,J ) =CITOLD  (ITR,J # ITP)  BR- 

180  CONTINUE  BR- 

WRITE(JW, 10009)  BR- 

KRITE ( JW ,10030)  PSENS  BR- 

IFdFINAL  .NE.  1)  WRITE  (JW  ,1 0029 ) MODE  BR- 

IFdFINAL  .EQ.  1)  WRITE  (JW  ,1 003  2)  MODEE  BR- 

IF (MODE  .EQ.  0 .OR.  I FINAL  .EQ.  1)  GO  TO  190  BR- 

IFTC-2  BR- 

IF(ITP  .EQ.  2)  IFTC-1  BR- 

WRITE(JW, 10031)  FTC (I  FTC)  BR- 

190  CONTINUE  BR- 

IFdFINAL  .NE.  1)  GO  TO  200  BR- 

PHIOBS-PHI1  BR- 

IF(ITP  .EQ.  2)  PHIOBS-PHI2  BR- 

WRITE(JW ,10035)  PHIOBS  BR- 

200  CONTINUE  BR- 

WRITEfJW, 10010)  BR- 

DEN-2.0*PI*RHO(1)  BR- 

IF(DEN.NE.0. .OR.MODE.NE.l) WRITE (JW ,10011)  RHO(l ), Z (1 ), CZERO, CZERO ,BR- 
SC  ZERO  BR- 

IF(DEN.EQ.0. .AND. MODE. EQ.l) WRITE (JW, 10012)  RHO(l) , Z (1 ), CZERO, ZER  ER- 

DO  210  1=1 , (NPTS-2)  BR- 

CALL  MAGPHS(CIT(I,J) ,CMP)  BR- 

DEN=  2. 0*P I*RHO( I+i)  BR- 

CITD=CIT(I,J)/CMPLX(DEN,0.0)  BR- 

A-CABS(CITD)  BR- 

210  WRITE(JW, 10013)  RHO(I+l),  Z(I+1),  CIT(I,J),  CMP,  CITD,  A BR- 

DEN=2.0‘PI*RHO(NPTS)  BR- 

IF(DEN.NE.0 . . OS, MODE. NE. i) WRITE (JW , 10011) RHO(NPTS) , Z(NPTS) , BR- 

5CZERO, CZERO, CZERO  BR- 

IF(DEN.EQ.0 , . AND. MODE. EQ.l) WRITE (JW,1 0012) RHO(NPTS) ,Z(NPTS) , CZERO, BR- 
SZER  BR- 

IFdFINAL  .NE.  1)  WRITE  (JW  ,1  00  21)  ACOND,  CD  ETRM , POWER  BR- 

IF(MODE  .NE.  0 .AND.  MlPR  .NE.  1)  WRITE (JW  ,1 0034)  ETl  BR- 

WRIT£(JW, 10014)  BR- 

WRITE(JW, 10030)  PSCNS  BR- 

IFdFINAL  .NE.  1)  WRITE  (JW  ,1 0029)  MODE  BR- 


4180 
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IF(IFINAL  .'EQ.  1)  WRITE (JW ,10032)  MODES  BR- 

IF{MODE  .EQ.  0 .OR.  I FINAL  .EQ.  1)  GO  TO  220  BR- 

IFTOl  BR- 

IF(ITP  .EQ.  2)  I FTC* 2 BR- 

WRITEdW  ,10031)  FTC  (I  FTC)  BR- 

220  CONTINUE  BR- 

IFdFINAL  .NE.  1)  GO  TO  230  BR- 

PHIOBS-PHI2  BR- 

I F ( I TP  .EQ.  2)  PHIOBS-PHI1  BR- 

WRITE(JW, 10035)  PHIOBS  BR- 

230  CONTINUE  BR- 

WRITE(JW, 10010)  ' BR- 

IS-NPTS-1  BR- 

IE«NPTS*2-3  BR- 

DO  240  I-IS,I£  BR- 

IR-I-1S+1  BR- 

CALL  MAGPHS(CIT(I,J) ,CMP)  BR- 

AMD- REAL (CMP)  BR- 

PH AS  E*  A I MAG (CMP)  BR- 

DEN-  2 . 0 *P I*  RHOPH ( I R)  BR- 

CITT*CIT(I ,J) *CMPLX(DEN,0.0 ) BR- 
A-CABS (CITT)  BR- 

240  WRITE(JW, 10013)  RHOPH (IR),  ZPH(IR) , CITT,  A,  PHASE,  CIT(I,J),  AMD  BR- 

IFdFINAL  .NE.  1)  WRITE  (JW  ,1 0021)  ACOND,CDETRM , POWER  BR- 

IF (MODE  .NE.  0 .AND.  MlPR  .NE.  1)  WRITE (JW ,1 0034)  EP1  BR- 

IF(IDB  .EQ.  0)  GO  TO  310  BR- 

WRITE(JW, 10015)  BR- 

WRITE(JW, 10030)  PSENS  BR- 

IF(IFINAL  .NE.  1)  WRITE(JW ,10029)  MODE  BR- 

IF (I FINAL  .EQ.  1)  WRITE (JW ,1 0032)  MODEE  BR- 

IF (MODE  .EQ.  0 .OR.  IFINAL  .EQ.  1)  GO  TO  250  BR- 

I FTC-1  BR- 

IF(ITP  .EQ.  2)  I FTC* 2 BR- 

WRITE(JW, 10031)  FTC (I  FTC)  BR- 

250  CONTINUE  BR- 

IF(IFINAL  .NE.  1)  GO  TO  260  BR- 

PHIOBS-PHI2  BR- 

IF(ITP  .EQ.  2)  PHIOBS-PHIi  BR- 

WRITE(JW, 10035)  PHIOBS  BR- 

260  CONTINUE  BR- 

WRITE(JW, 10010)  BR- 

DEN- 2. 0 *P I*  RHO(l ) BR- 

IF(DEN.NE.0 . .OR. MODE .NE. 1) WRITE/ JW ,10011)  RHO(l ) , Z (1 ) , CZERO,CZERO ,BR- 

$C  ZERO  BR- 

IF(DEN.EQ.0. .AND. MODE. EQ.l) WRITE (JW, 10012)  RHO(l) , Z(l) ,CZERO,ZER  BR- 
IS-IS+NPTS-1  BR- 
IE-1 E+NPTS-  2 BR- 

DO  270  I*1S,I E ER- 

IR-I-IS+2  BK- 

CALL  MAGPHS(CIT(I,J) ,CMP)  BR- 

DEN-2.0*PI*RHO(IR)  Bft- 

CITD*CIT( I ,J) /CMPLX(DEN,0 .0 ) BR- 
A-CABS (CITD)  BR- 

270  WRITE(JW, 10013)  RHO(IR) , Z(IR),  CIT(I,J),  CMP,  CITD,  A BR- 

DEN  »2.0*PI*RHO(NPTS)  BR- 

IF(DEN.NE.0. .OR,MQDE,NE,1)WRITE(JW,10011) RHO(NPTS) , Z (NPTS) , BR- 

$CZERO,CZERO,CZERO  BR- 

IF (DEN . EQ.0 . . AND, MODE.EQ. 1 ) WRITE (JW, 10012) RHO(NPTS) , Z(FPTS) ,CZER0,BR- 
SZ  ER  BR- 

IFdFINAL  .NE.  1)  WRITE  (JW ,1 0021)  ACOND,CDETRM, POWER  BR- 

IF(MODE  .NE.  0 .AND.  MlPR  .NE.  1)  WRITE (JW ,1 0034)  ET2  BR- 


4790 
4800 
4810 
4820 
4830 
4840 
4350 
4860 
4870 
4880 
4890 
4900 
4910 
4920 
4930 
494  0 
4950 
4960 
4970 
4980 
4990 
5000 
5010 
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518  0 

519  0 
5200 
5210 
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WRITE (JW, 10016)  BR- 

WRITE(JW, 10030)  PSENS  BR- 

IF(IFINAL  .NE.  1)  WRITE  (JW  ,10029)  MODE  BR- 

IF(IFINAL  ,EQ.  1)  WRITE  (JW  ,1 0032)  MOLEE  BR- 

IF(MODE  .EQ.  0 .OR.  IFINAL  .EQ.  1)  GO  TO  280  BR- 

IFT02  BR- 

IF(ITP  .EQ.  2)  I FTC* 1 BR- 

WRITE)JW ,10031)  FTC (I FTC)  BR- 

280  CONTINUE  l;R- 

IF{IFINAL  .NE.  1)  GO  TO  290  BR- 

PHI08S*PHI1  BR- 

IF (ITP  .EQ.  2)  PHIOBS-PHI2  BR- 

WRITE(JW ,10035)  PHIOBS  BR- 

290  CONTINUE  BR- 

WRITE(JW,  10010)  ’.'A- 

IS*  I S+NPTS-  2 E;<- 

DO  300  I=IS,NUNK  Bk»- 

IR-I-IS+1  r>  51- 

CALL  MAGPHS(CIT(I,J)  , CMP)  BR- 

AMD*REriL(CMP)  BR- 

PHASE“AIMAG (CMP)  BR- 

DEN*  2 . 0 * P I * RHOPH ( I R)  BR- 

CITT*CIT(I,J) *CMPLX{DEN,0.0)  BR- 

A*CABS (CITT)  BR- 

300  WRITE(JW, 10013)  RHOPH (IR) , ZPH(IR) , CITT,  A,  PHASE,  CIT(I,J),  AMD  BR- 

IFdFINAL  .NE.  1)  WRITE (JW ,1 0021)  ACOND,CDETRM , POWER  BR- 

IF(MODE  .NE.  0 .AND.  MlPR  .NE.  1)  WRITE(JW ,1 0034)  EP2  BR- 

310  CONTINUE  BR- 

320  CONTINUE  BR- 

IF(JM  .EQ.  NMODEM  .AND.  NMODEM  .NE.  1 .AND.  MODEB  .EQ.  0)  Bfc- 

51  FINAL  - IFINAL  + 1 BR- 

IFdFINAL  .EQ.  1)  GO  TO  120  BR- 

330  CONTINUE  BR- 

340  CONTINUE  BR- 

10000  FORMAT(///'0PO:NT'6X, 'RHO'laX, 'Z'llX, 'DELTA-T' 9X, 'GAMMA'//)  BR- 

10001  FORMAT( 14 , 1P2E15.6)  BR- 

10002  FORMAT (4 X, 1 P4E 1 5, 6 ) BR- 

10003  FORMAT( ' 1 ' 1 4 ( '*')/'  * '1 2X, ' * '/'  * INPUT  DATA  *'/'  *'12X,'*'/  BR- 

$'  '14( '*')///)  . BR- 

10004  FORMAT ( ' MAXIMUM  NUMBER  OF  MODES  TO  BE  USED  - '13//'  ORDER  OF  GAUSBR- 

5SIAN  QUADRATURE  * ’13/)  BR- 

10005  FGRMATC  FREE  SPACE  WAVELENGTH  - '1PE/)  BR- 

10006  FORMAT ( ' A PERFECT  CONDUCTOR  HAS  BEEN  ASSUMED.'/)  8R- 

10037  FORMAT) ' A DIELECTRIC  BODY  HAS  BEEN  ASSUMED: ' //6X, 'RELATIVE  DIELECBR- 

5TRIC  CONSTANT  » ' 1 PE//6  X, ' RELATIVE  PERMEABILITY  *>  '1PB/)  BR- 

10008  FORMATC  INCIDENT  FIELD  D ATA :' //6X ,' THETA  - ' 1PE//6X, 'PHI  * ' 1PE//BR- 

&6X, 'E-THETA  = '1 PE//6X, ' E-PHI  » ’1PE/)  BR- 

10009  FORMAT! 1 1 ' 4 3X, 1 "T"  COMPONENT  OF  ELECTRIC  CURRENT')  BR- 

10010  FORMAT ( 7 X,  'RHO'llX,  ' Z ' 1 0X, ' REAL' 7X, 'IMAGINARY' 4X, 'MAGNITUDE' 6X,  BR- 

$ 'PHASE' 8X, 'REAL'  7X, ' IMAGINARY' 4X, ' MAGNITUDE' /3 2X, 'TOTAL' 8X,  BR- 

$ 'TOTAL' 8X, 'TOTAL' 20X, 'DENSITY ' 6X, 'DENSITY ' 6X, 'DENSITY  ' /6X,  BR- 

$'  ***'10X, '***'9X, '****' 7X,9('*‘), 4 X, 9)'*'), 6 X, ******' 8X,  BR- 

$'•*** '7X,9 ('*') ,4X,9 ('*')//)  BR- 

10011  FORMAT) ' '1P5E13.5, ' 777  '1P3E13.5)  BR- 

10012  FORMAT) ' ' 1P5E13.5 ,1X,4 ( ' 777  '))  BR- 

10013  FORMAT) ' ’1P9E13.5)  . BR- 

10014  FORMAT) '1 '42X, ' "PHI"  COMPONENT  OF  ELECTRIC  CURRENT1)  BR- 

10015  FORMAT) '1 '43X, '"T"  COMPONENT  OF  MAGNETIC  CURRENT')  BR- 

10016  FORMAT) '1' 4 2X, '"PHI"  COMPONENT  OF  MAGNETIC  CURRENT')  BR- 

10017  FORMAT ( '0'17( '*')/'  *'15X, '*'/'  * COMPUTED  DATA  *'/'  * ’ 15X, ' * '/IX, BR- 

$17 {'*'>///)  BR- 


5400 
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10018 

FORM.*.?  (5 X,  ‘DIELECTRIC  CONSTANT  = ’1PE//5X,  'PERMEABII  ITY  = ' 

1PE//5XBR- 

6010 

WAVENUMBER  = '1 PE//5X, 'SPEED  OF  LIGHT  =*  ’1PE/) 

BR- 

6020 

10019 

FORMAT ('  PARAMETERS  OF  THE  DIELECTRIC  BODY:'/) 

BR- 

6030 

10020 

FORMAT ( ' OMEGA  = '1PE//'  FREQUENCY  = ’1PE/) 

BR- 

6040 

10021 

FORMAT(/'0INV.  COND.  NO. : ‘1  PE, 3X, ' DETERMINANT:  (* 

BR- 

6050 

$1P2E , ' ) TIMES  (10.)**('F6.1, ')*) 

BR- 

6060 

10022 

FORMAT ( 1 PARAMETERS  OF  FREE  SPACE:'/) 

BR- 

6070 

10023 

FORMAT ( 5 X, 'WAVELENGTH  = '1PE/) 

BR- 

6080 

10024 

FORMAT ( ' THE  MAXIMUM  NUMBER  OF  MODES  WAS  USED:’/’  THE  MAXIMUM 

(OVBR- 

6090 

$ER  NO.  OF  :NC.  FIELDS)  "RMS"  ERROR  IN  THE  CURRENT  VECTORS  WAS:' 

BR- 

6100 

S1PE/) 

BR- 

6110 

10025 

FORMAT ( ' A DELTA-GAP  SOURCE  IS  ASSUMED  AT  RHO  *'1PE,'  AND 

Z - ' 

BR- 

6120 

$1PE/) 

BR- 

6130 

10026 

FORMAT(46X, 'OUTPUT  IN  "MODE  1"  FORMAT'/) 

BR- 

6140 

10027 

FORMAT (8E) 

BR- 

6150 

10028 

FORMAT(8I) 

BR- 

bl6  0 

10029 

FORMAT ('0' 5 2X,' MODE  NUMBER  *12/) 

BR- 

6170 

10030 

FORM AT ( 3 9X, 'FOR  'A5,'  DIRECTED  PORTION  OF  INCIDENT  FIELD*) 

BR- 

6180 

10031 

FORMAT ( 4 9X, A3, ' (PHI)  COEFFICIENTS’/) 

FORMAT(  '0  ’4"7X,  'SUM  OF  MODES  0 THROUGH  *12/) 

BR- 

6190 

10032 

• BR- 

6200 

10033 

FORMAT ( '0WARNING:  TOO  MANY  INCIDENT  FIELDS  HAVE  BEEN  USED. 

ONLY  TBR- 

6210 

SHE  FIRST  11  HAVE  BEEN  PRINTED.') 

BR- 

6220 

18034 

FORMAT( ' 0ROOT-MEAN-SQUARE  ERROR  AT  THIS  STAGE:  '1PB) 

BR- 

6230 

10035 

FORMAT(45X, 'OBSERVED  AT  PHI  - 'F5.1,'  DEGREES'/) 

FORMAT ( '1'50X,27( '*')/51X, '*'25X,'*'/51X,' * CASE  NUMBER:' 

$ 73, 8X, ' * ' /5 IX, ' * *2  5X, ' * '/51X, ' * EXCITATION  NUMBER: *13,2 X, 

BR- 

6240 

10036 

BR- 

6250 

BR- 

6260 

S '*'/51X, '*'2  5X, ' * ' /51X,27 ( '*')/) 

BR- 

b2/0 

11100 

FORMAT  (1P2E,10X,1P2E) 

BR- 

6280 

11200 

FORMAT (3 ( ' CZ('I2I’,'I2,')»'1P2E,1X)) 

BR- 

6290 

RETURN 

BR- 

6300 

END 

bh- 

0310 

SUBROUTINE  ZMATRXiCZ , RHO, Z.RHOPH, ZPH,DBLTAT, GAMMA, NUNK ,M) 

BR- 

6320 

/•a*******************************"** ****************************** ******BR~ 

6330 

c 

BR- 

6340 

c 

SUBROUTINE  "ZMATRX"  FILLS  THE  IMPEDANCE  MATRIX  FOR  EACH  MODE. 

BR- 

6350 

c 

BR- 

6360 

L“J"L 

IMPLICIT  COMPLEX  (C) 

1 ”* “BR- 
BR- 

63/0 

6380 

COMPLEX  G1R,G2R,G4R,G5R 

BR- 

639  0 

COMPLEX  G0,G1,G2,G3,G5,G6 

BR- 

6400 

COMPLEX  G1M ,G  2M 

BR- 

6410 

REAL  MU1 , MU2 , COSG 

BR- 

6420 

DIMENSION  CZ (NUNK, NUNK) 

BR- 

6430 

DIMENSION  RHO (1 ) , Z (I ) , RHOPH (1 ) , ZPH (1 ) , DELTAT (1 ) , GAMMA (1 ) 

BR- 

6440 

CUMMON/PNT/NPTS,NUNKT,NUNKPH,NUNK2 

BR- 

6450 

common/int/nmodem,idb,ngq 

BR- 

6460 

COMMON/F  RQ/ P I , AK 1 , AK  2 , SL1 , SL2 , OMEG A 

BR- 

6470 

CUMMON/PRM/EPSl,EPS2 ,MUl,MU2 

BR- 

6480 

COMMON/FNC/RF , ZF,TF, RSL , ZSL , DEL , SING , COSG ,GM 

BR- 

6490 

COMMON/TP I/TWOP I , MAXP, ER, API 

BR- 

6500 

COMMON/S LF/ I SELF 

BR- 

6510 

EXTERNAL  G1 ,G 2 ,G 3 ,G 5 ,G 6 ,G0 ,G1R ,G 2R ,G4R ,G5R 

BR- 

6520 

EXTERNAL  G1M  ,G  2M 

BR- 

6530 

CHX(UUMMY) “CMP LX (DUMMY ,0.0) 

BR- 

6540 

CIX(DUMMY) “CMPLX(0.0, DUMMY) 

BR- 

6550 

NPTSM1 “NPTS-1 

BR- 

6560 

ER*l.E-2 

HR- 

6570 

NTE>NUHKT 

bh- 

0360 

MAXP=10000 

Bh- 

6390 

NPHEaNUNKPH 

BR- 

6600 

API=P I 

Bh- 

661  0 

9 7 


noon 


tm  FAOQt ISBSBT  mCIICABl* 

ib oi  oopy  ^Riusi^b  ttjmc 


TWOPI»2.0*PI  BR- 

XF(IDB  .EQ.  it)  GO  TO  10  BR- 

NTE*NTE/2  BR- 

NPHE-NPHE/2  BR- 

n'/'m^nte  br- 

NPHM*NPHE  BR- 

» 10  MODE=M  BR- 

> MODEP1-MODE+1  BR- 

nODEMl “MODE- 1 BR- 

NGQ2=NGQ*2  BR- 

KTS»NTE+NPHE  3R- 

KPHS*KTS+NTM  BR- 

Cti*«CMPLX(0.5,0.0)  BR- 

CZERO=>CNPLX(0. 0,0.0)  BR- 

CJO«CMPLX(0.0, OMEGA)  BR- 

CMU1=CMPLX( MU1,0 .0 ) BR- 

CMU2*CMPLX(MU2 ,0.0)  BR- 

CEPS1  =CMP1,X(  EPS1 ,0.0)  BR- 

CEPS2=CMPLX ( FPS2 , 0 ,0 ) BR- 

CO=*CMPLX  (OMPiGA,0 . 0 ) BR- 

CMP*CMrL..(  FLOAT(MODEPl)  ,0.0)  BR- 

CM*CMP1.X(  i’LOAT(MODE)  ,0.0)  BR- 

CMM^CMl  LX  i FLOAT  ( rlODEMl ) ,0.0)  BR- 

CM2“CM*CM  BR- 

BR- 

•EVALUATO  ALL  EXPRESSIONS  IN  WHICH  THE  FIELD  POINT  IS  LOCATED  AT  BR- 
A POSSIBLE  BEND  (NON-HALF  POINTS)  BR- 

BR- 

TL-0.0  AR- 

TK-0.5  BR- 

TU-1.0  LR- 

>J0  100  m»l,NTE  BR- 

IF3 «I Fl+KTS  BR- 

IF1P1-IF1+1  BR- 

RK-RKO(IFIPI)  BR- 

ZK*Z(IFlPl)  BR- 

DTK-DELTAT(IFl)  BR- 

DTKP1*DELTAT(IF1P1)  BR- 

sgk»sin(gamma;ifd ) br- 

SGKP1=S  I N (G  AMMA(  I F.lPl)  ) BR- 

ACGK=COS (G AMMA ( I Fl) ) BK- 

ACUKP1*C03 (GAMMA( IF1 Pi) ) BR- 

BRCKTS*  (DTKPl*SGKPl+DTK*SGK)/2,0  BR- 

BRCKTO (DTKPl*ACGKPl+DTK*ACGK)/2,0  BR- 

DO  90  JS1«1,NPHE  BK- 

152- IS1+NTE  BR- 

153- IS1+KTS  BR- 

IS4“I SI +KPHS  BR- 

ISlMi “I Sl-1  BR- 

IS3M1 “IS3-1  BR- 

RIM1 “RHO ( IS1)  BR- 

ZlMl-Z(ISl)  BR- 

t DTI= D E LT AT ( I S 1 ) BR- 

SGI»SIN(GAMMA(IS1) ) BR- 

ACGI“COS (GAMMA (I  SI) ) BR- 

RF-SK  BK- 

ZF»ZK  BK- 

RSL-RIM1  BR- 

ZSL=ZIM1  BK- 

DEL*DTI  BR" 

SING=SGI  BR- 


bb20 
bb30 
6640 
bb50 
bbb0 
bb  tit 
b680 
bb90 
b/00 
6710 
6720 
6730 
6740 
6750 
6760 
6770 
6780 
6790 
630  0 
6810 
6820 
6830 
6840 
6850 
6860 
6870 
6880 
6890 
690  0 
6910 
6920 
6936 
694  0 
6950 
6960 
6970 
6980 
6990 
7000 
7010 
/0  20 
7030 
7040 
/050 
/0b  V 
/0/0 
7080 
/B90 
7100 
/110 
7120 
.130 
/I  40 
7150 
/ 160 
U/0 
/lev 
/19  0 
uvv 
7210 
/ 220 
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nil  not  IS  BEET  QUAIrlT?  HttCWC ABLE 

nu  oopy  lumsupD  wipe 


COSG-ACGI 

C1B1 l“CIXtDTIwOMEGA*SGI*BRCKTS/2 .0 ) 
C3Bll=CIX(D'ri*OMEGA*ACGI*BRCKTC) 

C1B1 2*=CRX  ( -DTI*OMEGA*P I*BRCKTS) 

IF (IDB  .EQ.  0)  GO  TO  20 
C1B13*CIX(-DTI*ACGI*BRCKTS) 
C3B13-CIX(DTI*SGI*RK*BRCKTC) 
C5B13»CIX(-DTI*SGI*BRCKTS) 

C1B1 4*CRX ( DTl*TWOPI*BRCKTC) 
C2B14=CRX(-4.0*PI*DTI*RK*BRCKTC> 

C3B1 4=CRX( -DTl*TWOPI*BRCKTS) 

20  TF-TL 
ISELF*0 

IF( IS1  .EQ.  IFIP1)  ISELF=*1 
CALL  ELLPTC(TL,TM,NGQ,CANEL) 

CALL  ELLPTR(TL, TM,NGQ,CANELR) 

GM*FLOAT ( MODEPl ) 

CALL  CGQlT(G 1M  ,TL,TM ,NGQ,CANS , CANSR) 

CAlLP*CANS+CANt’,i. 

CAILPR’CANSR+CANELR 
GM=>  FLOAT  (MODEH1) 

CALL  CGQlT (GlM ,TL,TM,NGQ, CANS , CANSR) 

CA1LM-CANS+CANEL 
CA1 LMR=CANSR+CANELR 
IF(ISl  .EQ.  1)  GO  TO  30 
GM» FLOAT (MODE) 

CALL  CGQl (GlM ,TL,TM,NGQ, CANS) 

CA1  L=*CANS+CANEL 

CZ(IF1,IS1M1) =CZ(IF1,IS1M1) +ClBll*CMUl* (CA1LP+CA1LM) 
CZ  (I  FI , I S1M1 ) =CZ  ( IFl , I SI  Ml ) +C3BU*CMU1*CA1L 
IF(IDB  .EQ.  0)  GO  TO  50 

CZ (I F3 , I S3M1 ) =CZ (1 F3 , IS3M1 ) +C1B1 1*CEPS1 * (CA1 LP+CAl LM) 
C Z (I  F3 , 1 S3M1 ) *C Z (I  F3  , IS3M1 ) +C3B1 1*CEPS1 *CAlL 
30  I F ( I DB  .EQ.  0)  GO  TO  50 
GM- FLOAT (MODEPl) 

CALL  CGQlT (G2M ,TL,TM ,NGQ, CANS , CANSR) 

CA2LP=CANS+C ANEL 
CA2LPR=CANSR+CANELR 
GM»FLOAT(MODEMl ) 

CALL  CGQlT (G2M ,TL,TM,NGQ, CANS, CANSR) 

CA2LM=CANS+CANEL 

CA2LMR=CANSR+CANELR 

GM*FLOAT(MODE) 

IF( IS1  .EQ.  li  GO  TO  40 
CALL  CGQ1 (G  2M ;TL,TM  ,NGQ, CANS) 

CA2L-CANS+CANEL 

C Z (I  FI , I S1M1 ) =CZ ( I FI , IS1M1 ) +C1B1 1*CMU2* (CA2LP+CA2LM) 

CZ  (I  FI , IS1M1 ) =CZ(IF1,  IS1M1)  4C3BU*CMU2*CA2L 

CZ( I F3 , IS3M1 ) =CZ ( 1F3 , IS3M1 ) +C1B1 1*CEPS2* (CA2LP+CA2LM) 

CZ(IF3, IS3M1) =CZ ( IF3 , IS3M1 ) +C3Bll*CEPS2*CA2L 

CALL  CGQl (G0 , TL,TM ,NGQ,CANS0L) 

IF(1SELF  .EQ.  1)  CANS0L*CANS0L+CANALY(TL,TM,0,0> 

CALL  CGQl (G1,TL,TH  ,NGQ,CANS1L) 

CALL  CGQl (G3,TL,TM,NGQ, CANS3L) 

IF(ISELF  .EQ.  1)  CANS3L=CANS3L+CANALY(TL,TM,3,0) 

CB1 3=ClBl 3*CANS3L+C3B13*CANS0L+C5B13*CANS1L 
C 2 ( I FI , I S3M1 ) =C  Z (I Fl , I S3M1 ) + CB13 
40  CALL  CGQl (G4R ,TL,TM ,NGQ, CANS4L) 

CALL  CGQl (G5R  ,TL ,TM  ,NGQ, CANS 5L! 

IFdSELF  .EQ.  1)  CANS5L=CANS5L+CANALY(TL,TM,5,1) 

CALL  CGQ1(G2R,TL,TM,NGQ,CANS2L) 


BR-  /2  30 
BR-  7240 
BR-  7250 
BR-  7260 
BR-  7270 
BR-  7280 
BR-  729  0 
BR-  7300 
BR-  7310 
BR-  7320 
BR-  7330 
BR-  7340 
BR-  7350 
BR-  7360 
BR-  7370 
BR-  7380 
BR-  7390 
BR-  7400 
BR-  7410 
BR-  7420 
BR-  7430 
BR-  7440 
BR-  7450 
BR-  7460 
BR-  7470 
BR-  748  0 
BR-  7490 
BR-  /500 
BR-  7510 
BR-  7520 
BR-  7530 
BR-  7540 
BR-  7550 
BR-  7560 
BR-  7570 
BR-  7580 
BK-  /590 
BR-  7600 
BR-  7610 
BR-  7620 
BR-  7630 
BR-  7640 
BR-  /650 
BR-  7660 
BR-  7670 
BR-  7680 
BR-  7690 
BR-  770  0 
BR-  7710 
BR-  7720 
BR-  7730 
BR-  7740 
BR-  7750 
BR-  7760 
BR-  7770 
BR-  7780 
BR-  7790 
BR-  780  0 
BR-  /810 
BR-  7820 
BR-  7630 
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50  TF*TU  BR- 

ISELF-0  BR- 

IF(IS1  .EQ.  IF1)  ISELF*1  BR- 

CALL  ELLPTC(TM,TU,NGQ,CANEU)  BR- 

CALL  ELLPTR ( TM ,TU ,NGQ , CAN EUR)  8R- 

GM=*FLOAT(MODEPl)  BR- 

CALL  CGQ1T (G1M ,TM,TU ,NGQ,CANS ,CANSR)  BR- 

CA1UP*CANS+CANEU  BR- 

CA1  UPR*CANSR+CANF,UR  BR- 

GM=*FLOAT  ( MODEM1 ) BR- 

CALL  CGQ1T{G1M,TM,TU,NGQ,CANS,CANSR)  BR- 

CA1UM*CANS+CANEU  BR- 

CA1 UMR*CANSR+CANEUR  BRr 

IF(IS1  .EQ.  NPHE)  GO  TJ  60  BR- 

GM=FLOAT(MODE)  BR- 

CALL  CGQ1(G1M,TM,TU,NGQ,CANS)  BR- 

CA1UCANS+CANEU  BR- 

CZ (I Fl , I SJ.)  “CZ  (IFl , IS1)  +C1B1 1*CMU1*  (CA1UP+CA1 W)  +C3B11*CMU1*CA1U  BR- 
IF(IDB  .EQ.  0)  GO  TO  80  BR- 

CZ(IF3,IS3) =CZ(IF3, IS3) +C1B1 1*CEPS1 * (CAlUP+CAlUM) +C3Bll*CEPSl*CAlUBR- 
60  IF ( IDB  .EQ.  0)  GO  TO  80  BK- 

GM=FLOAT(MODEPl)  br- 

CALL  CGQ1T(G2M,TM,TU,NGQ,CANS,CANSR)  BR- 

CA2UP*CANS+C ANEU  BR- 

CA2UPR-CANSR+CANEUR  BR- 

GM»FLOAT(MODEMl)  BR- 

CALL  CGQ1T(G2M,TM,TU,NGQ,CANS,CANSR)  BR- 

C A2  UM=C ANS+C  ANEU  BR- 

CA2UMR-CANSR+CANEUR  BR- 

GM« FLOAT (MODE)  BR- 

IF( IS1  .EQ.  NPHE)  GO  TO  70  BR- 

CALL  CGQ1 (G2M ,TM  »TU ,NGQ ,CANS)  BR- 

CA2 U^CANS+CANEU  BR- 

CZ(IF1,  IS1)  »=CZ(IFl,  IS1)  +C1B1 1*CMU2*  (CA2UP+CA2UM)  +C3B1 1*CMU2*CA2U  BR~ 
CZ(IF3, IS3) =CZ(IF3,IS3) +C1B1 1*CEPS2* (CA2UP+CA2UM) +C3Bll*CEPS2*CA2UBR- 
CALL  CGQ1 (G0, TM ,TU ,NGQ, CANS0U)  BR- 

IFdSELF  .EQ.  1)  CANS0U«CANS0U+CANALY(TM,TU,0,0)  BR- 

CALL  CGQ1(G1,TM,TU,NGQ,CANS1U)  BR- 

CALL  CGQ1(G3,TM,TU,NGQ,CANS3U)  BR- 

IF( ISELF  .EQ.  1)  CANS3U-CANS3U+CANALY(TM,TU,3,0)  BR- 

CB1 3-ClBl 3"C ANS3U+C3B1 3*CANS0U+C5Bl 3*CANS1U  BR- 

C Z ( I Fl , I S3 ) =C  Z (I  Fl , I S3 ) + CB1  3 BR- 

70  CALL  CGQ1(G4R,TM,TU,NGQ,CANS4U)  BR- 

CALL  CGQ1(G5R,TM,TU,NGQ.CANS5U)  BR- 

IFdSELF  .EQ.  1)  CANS5U=*CANS5U+CANALY  (TM,TU  »5  >1 ) BR- 

CALL  CGQ1(G2R  ,TM  ,TU , NGQ, CANS 2U ) BR- 

CB14»C1B14*  (CANS4  L+C  ANS4U)  -‘■C2B1 4*  (CANS5L+CANS5U)  BR- 

CBl 4“CB1 4+C3B1 4* (CANS2L+CAN32U)  BR- 

C2(IF1, IS4) =CZ(IF1, IS4)  +CB14  BR- 

CA2MR*=CA2LMR+CA2UMR  BR- 

CA2PR=CA2LPR+CA2UPR  BR- 

CZdFl,  IS2)  =C Z (IF1 , IS2)  FC1B12«'CMU2*  (CA2PR-CA2KR)  BR- 

C Z (I F3 , 1 S4 ) =C  Z ( I F3 , 1 S4 ) +C1B1 2*CEPS2  * (CA2  PR-CA2MR)  BR- 

80  C1ftlMR=CAl  LMR+CAlUMR  BR- 

CA1  PR=*CA1  LPP.+CA1  UPR  BR- 

CZ(IF1 , IS2) =C  Z (IF1 , IS2)  +C1B1 2*CMU1* (CAlPR-CAlMR)  BR- 

IF(IDB  .EQ.  1)  CZ(IF3,  IS4) =CZ(IF3,  IS4 ) +C1B1 2*CEPS1* (CAlPR-CAlMR)  BR- 
90  CONTINUE  BR- 

100  CONTINUE  BR- 

BR- 

COMPUTE  ALL  EXPRESSIONS  IN  WHICH  THE  FIELD  POINT  IS  LOCATED  AT  BR- 


7840 
7850 
7860 
7870 
7880 
7890 
7900 
7910 
7920 
7930 
794  0 
7950 
7960 
7970 
7980 
/990 
8000 
8010 
8020 
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8040 
8050 
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80/0 
8080 
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8120 
8130 
8140 
8150 
8160 
8170 
8180 
819  0 
8200 
8210 
8220 
8230 
8240 
8250 
8260 
827  0 
8280 
8290 
8300 
8310 
83  20 
83  30 
8340 
8350 
8360 
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8380 
8390 
8400 
8410 
8420 
8430 
8440 
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A HALF  POINT  ON  THE  GENERATING  CURVE.  BR- 

Brm 

DO  270  IF1=»1,NPHE  BR- 

IF2-IF1+NTE  BR- 

IF3-IF1+KTS  BR- 

IF4*I Fl+KPHS  BR- 

IF1M1-IF1-1  BR- 

IF3Ml*IF3~l  BR- 

RK*RHOPH ( I FI ) BR- 

ZK-ZPH(IFl)  BR- 

DTK-DELTAT(lFl)  BR- 

SGK«SIN(GAMMA(IF1)  ) BR- 

ACGK=COS (GAMMA ( I FI ) ) BR- 

DO  260  IS1=1,NPHE  BK- 

IS2*IS1+NTE  BR- 

IS3“IS1+KTS  BR- 

IS4-IS1+KPHS  BK- 

IS1M1*IS1-1  BR- 

IS3m1=»IS3-1  BR- 

RlMl-RHO(ISl)  BR- 

ZIMi-Z(ISl)  BK- 

DTI’DELTAT(ISl)  BR- 

SGI“S IN (G AMMA( IS1)  ) BR- 

ACGI-COS(GAMMA(ISl)  ) BR- 

RF-RX  BR- 

ZF-ZK  BR- 

RSL-RIM1  BK- 

ZSL-Z1M1  BR- 

DEL-DTI  BR- 

SING-SG  BR- 

C0SG«»''JI  bk- 

C5B1I»CIX( 1.0/OMEGA)  BR- 

C2B1 2*CRX ( -DTI *TWOP I/OMEGA)  BR- 

C1B2  2*CIX( DTK*DTI*P I*OMEGA)  BR- 

C2B2  2«C I X( -DTK*DTI*TWOPI/ ( OMEGA *RK) ) BR- 

ClB21*CRX(DTK*DTI*OMEGA*SGl/2.0)  BR- 

C382 1*CRX ( -DTK/ (OMEGA*RK) ) BR- 

I F ( I DB  .EQ.  0)  GO  TO  110  BR- 

C1B23=CRX(-2.0»DTX*DTI*RK*ACGI)  BR- 

C3B2  3=»CRX  ( -DTK*DTI*ACGI)  BR- 

C5B23«CRX(DTK*DTI*SGI)  BR- 

C1B24«CIX(-DTK*TW0PI*DTI)  BR- 

110  TF-TM  BR- 

ISELF-0  BR- 

IF(IS1  .EQ.  IFl ) ISELF-1  BR- 

CALL  ELLPTC (TL, TM ,NGQ, CANEL)  BR- 

CALI  '•LLPT’P(TL,TK,NGQ,CANELR)  BR- 

Gha  Y L jAT ( MODEP 1 ) BR- 

CALL  CGQlT(GlM ,TL,TM,NGQ,CANS ,CANSR)  BR- 

CA1LP-CANS+CANEL  BR- 

CA1LPR=CANSR+CANELR  BR- 

GM=FL0AT(!.0DEM1)  Bft- 

CALL  CGQlV  (GlM  ,TL,TM  ,NGQ,  CANSf  CANSR)  AR- 

CA1LM=CANS^ANEL  BR- 

CAILMR=CANSR+CANELR  B.I- 

GM»FLOAT(MODE)  BR- 

CALL  CGQ1 (GlM ,TL,TM  ,NGQ,CANS)  BR- 

CA1L=CANS+CANEL  BR- 

IF(IS1  .EQ.  1)  GO  TO  120  BR- 

CZ(I?2,IS1M1) =CZ(IF2,ISlMl) +C1B2 1*CMU1* (CA1LP-CA1  IK)  BR- 

IF(IDB  .EQ.  0)  GO  TO  140  BR- 
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CZ (I F4 , IS3M1 ) =CZ ( IF4  f I S3M1 ) +C1B2 l*CEPSl * (CA1LP-CA1LM) 
120  IFUDB  .EQ.  0)  GO  TO  140 
GM“FLOAT  ( MODEPl ) 

CALL  CGQ1T (G 2M , TL , TM ,NGQ , CANS , CANS R) 

CA2LP“CANS+CANEL 
CA2LPR=CANSR+CANELR 
GM“FLOAT ( MODEM1 ) 

CALL  CGQ1T (G2M ,TL, TM,NGQ, CANS ,CANSR) 

C A2  LM=C  ANS+  C AN  EL 

CA2LMR=*CANSR+CANELR 

GM=FLOAT(MODE) 

CALL  CGQl(G2M ,TL,TM  ,NGQ,CANS) 

CA2L-CANS+CANEL 
IF(IS1  .EQ.  1)  GO  TO  130 

CZ(IF2, IS1M1) =CZ(IF2,IS1M1) +C1B21*CM02* (CA2LP-CA2LM) 
CZ(IF4,IS3M1)  =CZ(IF4f  IS3M1 ) +C1B21*CEPS2*  (CA2LP-CA2LM) 
CALL  CGQ1(G5,TL,TM ,NGQ,CANS5L) 

IFdSELF  .EQ.  1)  CANS5L«CANS5L+CANALY <TL,TM,5 , 0 ) 

CALL  CGQ1(G6,TL,TM,NGQ,CANS6L) 

CALL  CGQ1 (G2,TL,TM ,NGQ, CANS2L) 
CB23“C1B23*CANS5L«MB23*CANS6L+C5B23*CANS2L 
CZ ( I F2 , IS3M1 ) =CZ( IF2 , IS3M1 ) + CB23 
130  CALL  CGQ1 (G1R ,TL,TM ,NGQ, CANS1L) 

140  CONTINUE 

CALL  ELLPTC(TM,TU,NGQ,CANEU) 

CALL  ELLPTR(TM,TU,NGQ,CANEUR) 

GM- FLOAT (MODEPl) 

CALL  CGQ1T (GlM ,TM ,TU ,NGQ , CANS , CANS R) 

CAlUP-CANS+CANEU 

caiupr*cansr+caneur 

GM-FLOAT(MODEMl) 

CALL  CGQ1T(G1m ,TM ,TU ,NGQ, CANS , CANSR) 

CA1  UM=*CANS+C  ANEU 
CA1 UMR=CANSR+CANEUR 
GM” FLOAT (MODE) 

CALL  CGQ1 (GlM ,TM,T0,NGQf CANS) 

CAl U-CANS+CANEU 
CA1.-CA1L+CA1U 
CA1PR»CA1LPR+CA1UPR 
CA1MP»CA1 LMR+CA1UMR 
IFdSl  .EQ.  NPHE)  GO  TO  150 

C 2 (I F2  , IS1)  “CZ  ( I F2 , ISl)  +C1B21*CMU1*  (CA1UP-CA1UM) 
IF(IDB  .EQ.  0)  GO  TO  230 

CZ  1 1 F4 , 1S3 ) “CZ(IF4,  IS3)  +C1B21*CEPS1*  (CA1UP-CA1UM) 

150  IF( IDB  , EQ.  0)  GO  TO  230 
GM»FLOAT( MODEPl) 

CALL  CGQ1T (G2M ,TM ,TU ,NGQ, CANS ,CANSR) 
ca2up“Cans+ca.;eu 
CA2  UPR“CANSR+CANEUR 
GM“FLOAT(MODEMl ) 

CALL  CGQ1T (G2K ,TM,TU ,NGQ , CANSjCANSR) 

CA2UM-CANS+CANEC 
CA2UMR=CANSR+CANEUR 
GM» FLOAT (MODE) 

CALL  CGQ1 (G2M ,TM ,TU,NGQ, CANS) 

CA2U“CANS+CANEU 

IF(IS1  .EQ.  NPHE)  GO  TO  160 

CZ(IF2,IS1) =CZ ( IF2 , IS1) +C1B21*CMU2* (CA2UP-CA2UH) 

CZ (IF4 / IS3) =CZ ( IF4 , 1S3) +C1B2 1*CEPS2* (CA2UP-CA2UM) 
CALL  CCQ1 (G5 ,TM,TU#NGQ, CANS5U) 

IF(ISELF  .EQ.  1)  CANS5U=CANS5U+CANALY(TM,TU,5,0) 
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BR-  9500 
BR-  9510 
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CALL  CGQ1 (G6,TM,TU,NGQ,CANS6U) 

CALL  CGQ1<G2,TM,TU,NGQ,CANS2U) 
CB23=C1B23*CANS5U+C3B23*CANS6U+C5B23*CANS20 

CZ(IF2,IS3) -CZ(IF2,IS3)  + CB23 
160  CALL  CGQ1(G1R,TM,TU,NGQ,CANS1U) 

CZ(IF2, IS4) =CZ {IF2, IS4) +C1B24* (CANS1L+CANS1U) 

CA2=CA2L+CA2'J 

CA2PR*CA2LPR+CA2UPR 

CA2MR=CA2LMR+CA2UMR 

CA2EPS=CA2/CEPS2" 

CAMU=CAl/CMUl+CA2/CMU2 

CZ (IF2 , IS2) “CZ (IF2, IS2)  +C1B22*CM02* (CA2PR+CA2MR) 

CZ (IF2 , IS2 1 “CZ (IF2 , IS2) +C2B22*CA2EPS*CM2 
CB4  4=C1B2  2*  (CEPSl*  {CAl  PR+CA1MR)  +CEPS2*  (CA2 PR+CA2MR) ) 
CB44=CB44+CM2*C2B22*CAMU 
CZ(XF4,IS4)«CZ(IF4,IS4)  + CB44 
CBl 2“C2B1 2*CM*CA2EPS 
CB3  4=C2B1 2*CM*CAMU 
CB21»C3B21*CM*CA2EPS 
CB4  3XC3B2 1*CM*CAMU 
CBl lxC5Bl 1*CA2EPS 
CB33«C5Bll*CAMU 
IF  (I  FI  .EQ.  1)  GO  TO  170 
CZ(IF1M1,IS2)*CZ(IF1M1,IS2)  + CB12 
CZ ( IF3M1 , IS4) “C  Z ( IF3M1 , IS4)  + CB34 
170  I F ( I FI  .EQ.  NPHE)  GO  TO  180 

CZ (IF1 , IS2) “CZ(I Fl , IS2)  - CB12 
CZ(IF3,IS4) -CZ(IF3,IS4)  - CB34 
180  IFdSl  .EQ.  1)  GO  TO  200 

CZdF2,ISlMl) -CZ(IF2,iSlMl)  - CB21 
CZ ( IF4  , IS3M1 ) “CZ  ( I F4  , IS3M1 ) -CB43 
IF ( I FI  .EQ.  1)  GO  TO  190 
C Z ( I F1M1 , IS1M1 ) “CZ ( I FlMl » ISlMl ) ~ CB11 
CZ (IF3M1 , IS3M1) =CZ (IF3M1 * IS3M1)  - CB33 
190  IF (I FI  .EQ.  NPHE)  GO  TO  200 

CZ(IF1,IS1M1)«CZ(IF1,IS1M1)  + CB11 
CZ ( IF3 , IS3M1) “CZ ( IF3 , IS3M1)  + CB33 
200  IF(IS1  .EQ.  NPHE)  GO  TO  220 

CZ{IF2,IS1)»CZ<IF2,IS1)  + CB21 
CZ(IF4,IS3)-CZ(IF4,IS3)  + CB43 
IF(IF1  .EQ.  1)  GO  TO  210 
CZ(IF1M1,IS1}*CZ(IF1M1,IS1)  + CB11 
CZ (IF3H1 , IS3)  “CZ  (IF3M1 » IS3)  +CB33 
210  I F ( I Fl  .EQ.  NPHE)  GO  TO  220 
CZ(IF1,IS1)«CZ(IF1,IS1)  - CB11 
CZ(IF3  , IS3)  «CZ  (IF3  ( IS3)  - CB33 
220  CONTINUE 
230  CONTINUE 

CA1EPS=CA1/CEPS1 
CB21»C3B21*CM*CA1EPS 
CB12=C2Bl 2*CM*CA1EPS 

CBl l“C5Bl 1* CAl EPS  . , 

CZ(IF2,IS2)*CZ(IF2f IS2) +C1B22*CMU1* (CA1PR+CA1RR) 

CZ(IF2,IS2)=CZ(IF2,IS2)+C2B22*CA1EPS*CM2 

IF (IF1  .NE.  1)  CZ(IF1M1, IS2) *CZ (I F1M1 , IS2)  + CB12 
X F (I Fl  .NE.  NPHE)  CZ d Fl , IS2) «CZ d Fl,  IS2)  - CB12 
IFdSl  .EQ.  1)  GO  TO  240 
CZ  dF2 / IS1M1 ) =CZ  d F2  , IS1M1 ) - CB2  1 

IF(IF1  .NE.  1)  CZ (IF1M1 , IS1N1 ) *CZ (IF1M1 » IS1M1 ) - CBil 
I F ( I Fl  .NE.  NPHE)  CZ  (I  Fl , IS1M1 ) =>CZ  (I  Fl , IS1H1 ) + CB11 

240  IFdSl  .EQ.  NPHE)  GC  TO  250 


Bit*  9670 
BR-  9680 
BR-  9690 
BR-  9700 
BR-  9/10 
BR-  9720 
BR-  9730 
BR-  9740 
BR-  9750 
BR-  9760 
BR-  9770 
BR-  9780 
BR-  9790 
BR-  9800 
BR-  9810 
BR-  9820 
BR-  9830 
BR-  984  0 
BR-  9850 
BR-  9860 
BR-  9870 
BR-  9880 
BR-  9890 
BR-  9900 
BR-  9910 
BR-  9920 
BR-  9930 
BR-  9940 
BR-  9950 
BR-  9960 
BR-  9970 
BK-  9980 
BR-  9990 
BR-10000 
BR- 10010 
BR-10020 
BR-10030 
BR-10040 
BR-10050 
BR-10060 
BR-10070 
BR-10080 
8R-10090 
BR-10100 
BR-10110 
BR-10120 
BR-10130 
BR-10140 
BR-10150 
BR-1016  0 
BR-10170 
BK-10180 
BR-10190 
BR-10200 
BR-10210 
BR-10220 
BK-10ZJ0 
BR-10240 
BR-10250 
BK-1826  0 
BR-1027  0 


103 


rroroo  rnonoono 


IBIS  PAGE  IS  BEST  QUALITY  FRACTICABU| 
' Inoi  ODrXfUW^SHED  TO  DOC 


CZ ( IF2 , IS1) =CZ (IF2, IS1)  + CB21 

I F { X FI  ,NE.  1)  CZ ( XFlMl / IS1)  =CZ (I F1M1 , IS1)  + CB11 
IF(IF1  .NE.  NPHE)  CZ <IF1 , 181) -CZ (IF1 , IS1)  - CB11 
250  CONTINUE 
260  CONTINUE 
270  CONTINUE 


•NOW  WE  HAVE: 


(BETA- 31) 
(BETA- 32) 
(BETA- 41) 
(BETA- 42) 


» -(BETA-13) 
= -(BETA-14) 
* -(BETA-23) 
= -(BETA- 24) 


IF(IDB  .EQ.  0)  GO  TO  290 

NE»NTE+NPHE 

NE1=NE+1 

DO  280  IF*1,NE 

IFP=T F+NE 

DO  280  IS*NE1,NUNK 
ISP*! S-N  E 

CZ ( I FP,  ISP)  =*  -C  Z ( I F , I S ) 
280  CONTINUE 
290  CONTINUE 
RETURN 
EnD 
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SUBROUTINE  CVFILL(CVTHC ,CV?HC ,RHO , Z,RROPH, ZPH,DELTAT, GAMMA,  BR-10540 

# THETA , PHI , ETHETA, EPH I , ANTFD , NUNK ,NFLDS , MODE ) BR-10550 

***************************  ***W *»*»*»*»  ***  ************  *******  ******  **  *BR“  10560 


SUBROUTINE  CVFILL  COMPUTES  THE  FORCING  FUNCTION  VECTOR  (OR 
MATRIX  FOR  MULTIPLE  EXCITATIONS) . 


IMPLICIT  COMPLEX  (C) 

DOUBLE  PRECISION  BARG1 , BJ (50) , BY (2) 

REAL  COSTH,CGK ,CGKPl , MU 

DIMENSION  CVTHC (NUNK,NFLDS)  ,CVPHC(NUNK,NFI,D6) 

DIMENSION  RHO(l) ,Z(1) ,RHOPH(l) , ZPB(l) ,GAMMA(1) ,DELTAT{1) 

DIMENSION  THETA(l) ,PHI(1) , ETHETA (1) ,EPHI(1) , ANTFD (1) 

COMMON/PNT/NPTS , NUNXT , N UNK  PH, N UNK  2 

COMMON/ I NT/N MODEM , IDB ,NGQ 

COMMON/FRQ/PI , AK1, AK2, SLl , SL2 ,OMBGA 

CALL  2ERO(CVTHC,NUNK«NFLDS) 

CALL  ZERO(CVPHC,NUNK*NFLDS) 

MODEP1-MODE+1 

HODEM1 «MODE-l 

NTE»NUNKT 

NPHE-NUNKPH 

NTM*0 

NPHM»0 

IF(IDB  .EQ.  0)  GO  TO  10 
NTE-NTE/2 
NPHE*NPHE/2 
NTM*NTE 
NPHM=NPHE 
10  CONTINUE 
TPI»2.0*PI 
FPI*4.0*PI 
JMPI*MODEPl+l 
JM*MODEPl 


BR-10570 
BR-10580 
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BR-10600 
BR-10610 
BR-10620 
BR-1063  0 
BR-1064 • 
BR-10650 
BR-10668 
BR-10670 
BR-10680 
BR-10690 
BR-10700 
BR-10/10 
BR-10720 
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BR-10740 
BR-10750 
BR-10760 
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BR-1083  0 
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JMM1-M0DE 
MU-FPI‘l.E-7 
ETA-0  MEG A*MU/AK 1 
CJ«CMPLX(0. 0,1.0) 

EVALUATE  FIEI^S  AT  NON-HAL?  POINTS 

DO  110  IF1-1,NTE 
IF3  -IF1+NTE+NPHE 
IF1P1-IF1+1 
RF-HHO(lFlPl) 

ZF-Z(IFlPl) 

DTK-DELTAT(IFl)  . 

DTKPl-DELTAT(IFlPl) 

SGK-S IN {GAMMA (IF1) ) 

SGKP1«SIN(GAMMA(IF1P1)) 

CGK-COS (GAMMA ( I FI ) > 

CGKP1«C0S(GAMMA(IF1P1)  ) 

BRCKTS- (DTKPl*SGKPl+DTK*SGK)  /2.0 
BRCKTC- (DTKP1*CGK?1+DTK*CGK) ^2.0 
DO  100  I«),NFLDS 
COSTH-COS (THETA ( I ) ‘PI/180.0) 

SINTH-SIN(THETA(I) ‘PI/180. 0) 

BAP.G-AKl*RF*SINTH 
BARG1-DBLE (BARG) 

IF (BARG)  20,  20,  40 
20  DO  36  K»2,JMP1 
30  BJ (K) *0.0 
Bail) -1.0 
GO  TO  50 

40  CALL  BESSEL (BARGl,MODEPl ,BJ ,BY,1) 

50  BETl-TPI*COSTH*BRCKTS* ETHETA(X) 

BET2— FPI*SINTH‘BRCKTC*ETHETA(I) 
BET3-TPI*BRCKTS‘EPHI(I) 

BHTl-TPI*BRCKTS*COSTH*  EPHI (I) /ETA 
BHT2— FPI‘BRCKTC*SINTH*EPHI(I)  /ETA 
BHT3— TPI‘BRCKTS*ETHETA(I)/ETA 
CJMP1“CMPLX( SNGL (BJ (JMPl) ) ,1.0) 
CJM»CMPLX(SNGL(BJ(JM) },0.0) 

IF(JMMl)  60,  60,  70 
60  CJMM1 "-CJMP1 
GO  TO  80 

70  CJMM1«CMPLX(SNGL(BJ(JNM1) ) ,0.0) 

80  CONTINUE 

CJPMMi-CJ”  (MODEM1-1) 

CJPMMl*CJ*CJPMM2 

CJ?M*CJ*CJPMM1 

CJPMP1»CJ«CJPM 

EARG-AKl*ZF‘COSTH-PHI( I) ‘FLOAT(HODE) *P 1/180.0 
CARU-CMPLX(Z ,0 , EARG) 

CE=CEXP(CARG) 

CV*  (CJMM1*CJPMM1  +CJPMP1*CJMP1)  ‘CMPLX  (BET1 ,0.0) 
CV“CV+CJPM*CJM»CMPLX(BET2,0.0) 

CVTHCUFl  ,1)  =*CV*CE 

CV*  (CJPMM2*CJMM1-CJPM‘CJMP1)  ‘CMPLX  (BBT3  ,0,0  ) 
CVPHC ( I FI , 1) -CV*CE 
IF(IDB)  90,  100,  90 

90  CV»(CJPMM2*CJMKi-CJPM‘CJMPl)  ‘CMPLX  (BBT3., 0.0) 
CVTHC (I F3 , I) “CV"CE 

CV=* (CJ  PMM1*CJKK1  +CJPMP1*CJMP1) ‘CMPLX (BHT1 ,0.0) 
CV*CV+CJPM*CJM‘CMPLX(BHT2,0,0) 
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CVPHC(IF3,I) =CV*CE 

BR-11500 

100 

CONTINUE 

BR-11510 

110 

CONTINUE 

BR-11520 

BR-11530 

■ EVALUATE  FIELDS  AT  HALF  POINTS 

BR-11540 

BR-11550 

DO  200  IF1*1 ,NPHE 

BR-11560 

IF2*I Fl+NTE 

BK-115/0 

IF4  =1 F2  +NPHE+NTM 

BR-11580 

RF*RHOPH(IFl) 

BR-11590 

ZF*ZPH(IF1) 

BR- 11600 

DTK*DELTAT(IF1) 

BR-11610 

DO  190  I=1,NFLDS 

BR-1162  0 

COSTH=COS(THETA(I) *PI/180.0) 

BR-11630 

SINTH*SIN(THETA(I) *PI/180.0) 

BR-11640 

BhRG*AK1*RF*SINTH 

BR-11650 

BARGl=DBLE(BARGj 

BR-11660 

IF(BARG)  120,  120,  140 

BR-11670 

120 

DO  130  K«2,JMPl 

BR-11680 

130 

BJ (K) *0.0 

BR-11690 

BJ(1) *1.0 

BR-11700 

GO  TO  150 

BR-11710 

140 

CALL  BESSEL (BARG1 ,M0DEP1 ,BJ , BY, 1) 

BR-11720 

150 

BEP1»-DTK*TPI*C0STH*ETHETA(I) 

BR-11730 

BEP2=DTK*TPI*EPHI(I) 

BR-117  40 

BHP1—  DTK*TPI*COSTH*EPHI  (I)  /ETA 

BR-11750 

BHP2*-DTK*TPI*ETHETA{I)/ETA 

BR-11760 

CJMPl*CMPLX{ SNGL(BJ ( JMP1) ) ,0.0) 

BR-11770 

IF(JMMl)  160,  160,  170 

BR-117 80 

160 

CJMMl—CJMPl 

BR-11/90 

GO  TO  180 

BR-11800 

170 

CJMRl *CMPLX( SNGL(BJ ( JMMl) ) ,0.0)  ‘ 

BR-11810 

180 

CONTINUE 

BR-118  20 

CJPMM2 «CJ*  * (MODEM1-1 ) 

BR-11830 

CJPMM1 =CJ*CJPMM2 

BR-11840 

CJPM*CJ*CJPMM1 

BR-118  50 

CJPMPl=CJ  PM*CJ 

BR-11860 

EARG=AK1*ZF*C0STH-PHI(I) *FLOAT(MODB) *P 1/180. 

BR-118 70 

CARG*CMPLX(0 .8 , EARG) 

BR-11880 

CE-CEXP(CARG) 

BR-11890 

CV*  (CJPMM2*CJMM1-CJPM*CJMP1)  *CMPU  (BEP1,0.») 

BR-11900 

CVTHC (I F2 , I) »CV*CE 

BR-11910 

CV* (CJPMM1*CJMM1 +CJPMPl*CJMPl) *CMPLX (BEP2,0 .0 ) 

BR-119  20 

CVPHC ( I F2 , I) *CV*CE 

BR-11930 

IF(IDB  .EQ.  0)  GO  TO  190 

BR-119  40 

CV*  (CJPMM2*CJMM1-CJPM*CJMP1)  *CMPUC  (BHPl  ,0 .0 ) 

BR-11950 

CVPHC (I F4 , I) *CV*CE 

BR-119 60 

CV* (CJPMM1*CJMM1 +CJPMPl*CJMPl) *CMPLX (BHP2 ,0 .0 ) 

BR-11970 

CVTHC {I F4 , 1) *CV*  C E 

BR-119 80 

190 

CONTINUE 

BR-11990 

200 

CONTINUE 

BR-12000 

BR-12010 

- INCLUDE  VOLTAGE  SOURCE  IF  PRESENT 

BR-12020 

BR-12030 

IF(MODE  ,NE.  0)  GO  TO  220 

BR-12040 

DO  210  I* 1 , NFLDS 

BR-12050 

IANTFD=I FIX (ANTFD (I)  ) 

BR-12060 

IFdANTFD  .LE.  0)  GO  TO  210 

BR-12070 

CVTHC ( I ANTFD, I ) =CVTHC ( I ANTFD , I ) +FPI 

BR-12080 

210 

CONTINUE 

BR- 12090 

2 20 

CONTINUE 

BR-1213  0 
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RETORN 

END 

SUBROUTINE  SUMOUT(CI ,NUNK , NFL DS , MODE , I THETA) 


SUBROUTINE  “SUMODf  TRANSFORMS  THE  EXPONENTIAL  COMPONENT  CURRENT 
COEFFICIENTS  TO  SINE  AND  COSINE  COMPONENT  CURRENT  COEFFICIENTS. 


IMPLICIT  COMPLEX  (C) 

DIMENSION  Cl  (NUNK.NFT.DE) 

COMMON/ INT/NMODEM ,1  LB ,NGQ 
COMMON/PNT/NPTS,NUNXT,NUNKPH,NONK2 
IF (MODE  .EQ.  0)  RETURN 
NTE*N UNKT 

nphe-nunkph 

NTM-0 
NPHM-0 

IF (IDB  .EQ.  0)  GO  TO  II 
NTE-NTE/2 
NPHE^NPHE/2 
NTM-NTE 
NPHM“NPHE 
10  CONTINUE 

CF1 -CMP LX(2. 0,9.0} 

CF2-CMPLX(0. 0,2.0) 

IFIITHETA  .EQ.  1)  GO  TO  20 
CF-CP1 
CF1-CF2 
CF2-CF 
20  CONTINUE 

DO  60  J-l ,NFLD6 
DO  50  I«1,NUNX 
IF(I  .GT.  NTE)  GO  TO  41 
30  CI(I,J)-CI(I,J) *CPl 
GO  TO  50 

40  IF(I  .GT.  NTE+nPHE+NTM)  GO  TO  31 
CIU,J)-CI(I,J)  *CF2 
50  CONTINUE 
60  CONTINUE 
RETURN 
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END  BR-12320 

Subroutine  rmserr(Cit,citold,nunk,rmshxb)  br-12531 

»*****»*«»  »* *«••**•****•****»**••*•***•••• »*»»**»****»****»****»»»*  w*  *gp . 12540 


SUBROUTINE  ” RMSERR"  COMPUTES  THE  ROOT-MEAN-SQUARE  RELATIVE 
ERROR  BETWEEN  THE  PRESENT  CURRENT  SOLUTION  AND  THE  PREVIOUS 
SOLUTION. 


•no************************** *» *»*»** *****»»**?•* ****** ******** 

IMPLICIT  COMPLEX  (C) 

DIMENSION  CIT(NUNK)  , CITOLD(NUNK) 

SQERR-0.0 
DO  10  1*1, NUNK 
DEN-CABS (CIT(  I)  ) 

ANUM-CABS (CIT{ I) -CITOLD(I)  ) 

IF(DEN  .LT.  l.E-9)  DEN’CABS(C ITOLD (I ) ) 

IF(DEN  .LT,  l.E-9)  DEN-1.0 
RELERR-ANUM/DEN 
S0ERR=SQERR+RELBHR**2 
10  CONTINUE 


BH-12550 
.iR-12560 
BR-12570 
BR-12580 
BR-12300 
BR- 12600 
BR-12610 
BR-12620 
BR- 12630 
BR-12640 
BR-12650 
ER-12660 
BR-12670 
BR- 12680 
BR-12690 
BR-127  00 
BR-12710 


107 


nnnooor:  noon  o nr 


HSSSS&’“r* 


SQMERR»SQ£RR/FLOAT (NUNK) 

RMSMXE=SQRT(SQMERR) 

RETURN 

END 

SUBROUTINE  ZERO(C,N) 

INITIALIZATION  ROUTINE 

COMPLEX  C (N) 

C0 -CMP LX{0. 0,8.0) 

DO  10  J-1,N 
10  C(J)=*C0 
RETURN 
END 

SUBROUTINE  MAGPHS (CRI , CMP) 

ROUTINE  TO  CALCULATE  MAGNITUDE  AND  PHASE  OP  A COMPLEX  NUMBER 

RESULT  IS  STORED  IN  A COMPLEX  NUMBER  AS  HELL 

( 

IMPLICIT  COMPLEX  (C) 

REAL  I M 

COMMON/FRQ/PI,A,B,D,E,P 

R-REAL(CRl) 

I-AIMAG(CRI) 

M-CABS (CRI) 

P-ATAN2 ( I , R) *1 80.0/PI 
CMP«CMPLX(M,P) 

RETURN 

END 

SUBROUTINE  NORMAL (C I , ETHETA ,BPHI ,NUNK ,NPLD6 ) 


SUBROUTINE  NORMAL  NORMALIZES  THE  ELECTRIC  CURRENT  BY  THE 
INCIDENT  H-FIELD  AND  THE  MAGNETIC  CURRENT  BY  THE  INCIDENT 
E-FIELD. 


IMPLICIT  COMPLEX  (C) 

DIMENSION  01 (NUNK,NFLDS) 

DIMENSION  ETHETA (1) ,EPHI(1) 

COMMON/F  RQ/P I , AK 1 , AX  2 , SLl , SL2 , OMEGA 

COMMON/INT/NMODEM,IDB,NGQ 

NUNKE-NUNK 

IF(IDB  .EQ.  1)  nUNKE-NUNK/2 

NUNKE1-NUNKE+1 

DO  30  K-l , NFLDS 

ENORM-SQRT( ETHETA (K) **2+EPHI(K) **2) 
HNORM“ENORM*SLl*  8 ,8  5418*13  3E-12 
CEN»CMPL,X(ENORM,0.0) 

CHN-CMPLX (HNORM ,0.0) 

DO  10  I* 1 , NUNKE 
10  CI(I,K}=CI(I,K)/CHN 

IF(IDB  .EQ.  0)  GO  TO  30 
DO  20  I-NUNKEl.NUNK 
20  CI(I,K)-CI(1,K)/CEN 
30  CONTINUE 
RETURN 
END 

COMPLEX  FUNCTION  GIM(TP) 

IMPLICIT  COMPLEX  (C) 

COMMON/TP I/TKOP I , MAXP , ER, PI 
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COMMON/ARC/T 
EXTERNAL  CG1M 
T“TP 

CALL  TRPADP ( CG 1M  , 0 . 0 , P I , ER, MAX  P , I ER , C ANS) 

GlM-C ANS/C  MP  LX ( TWOPI ,0.0) 

RETURN 

END 

COMPLEX  FUNCTION  CGIM(PUIP) 

IMPLICIT  COMPLEX  (C) 

REAL  COSG 
COMMON/ARC/TP 

COMMON/FNC/RF , ZF ,TF, RSL , ZSL,DEL , SING , COSG ,GM 
COMMON/FRQ/P I , AK 1 , AK  2 , SL1 , SL2 , OMEGA 
RS- RSL+TP* DEL*  SING 
ZS“ZSL+TP* DEL* COSG 

R“SQRT(RF*RF+RS*RS+ (ZF-ZS) **2-2.0*RF*RS*COS(PHIP) ) 
CARG«CMPLX(0 . 0 ,-AKl*R) 

CR«CMPLX(R,0.0) 

ACMP-COS(GM*PHIP) 

CGiM“ (CMPLX( ACMP,0 . 0 ) *CEXP (CARG) -CMPLX(1. 0,0.0) )/CR 

RETURN 

END 

COMPLEX  FUNCTION  G2M(TP) 

IMPLICIT  CuMFLEX  (C) 

COMMON/T  P I/T WOP I , MAX  P , ER, P I 

COMMON/ARC/T 

EXTERNAL  CG2M 

T«TP 

CALL  TRPADP (CG2M ,0.0 ,P I ,ER, MAXP, IER,CANS) 

G 2M»C  ANS/C  MP  LX • TWOP 1 , 0 . 6 ) 

RETURN 

END 

COMPLEX  FUNCTION  CG2M(PHIP) 

IMPLICIT  COMPLEX  (C) 

REAL  COSG 
COMMON/ARC/TP 

COMMON/FNC/RF, ZF,TF,RSL,ZSL, DEL, SING, COSG,GM 
COMKON/FRQ/PI , AK1 ,AK2,SL1 , SL2 , OMEGA 
RS-RSL+TP*DEL*SING 
ZS»ZSL+TP*DEL*COSG 

R*SQRT(RF*RF+RS*RS+  (ZF^-ZS)  **2-2.e*RP*RS*COS(PHIP)  ) 
CARG-CMPLX(0.0,-AK2*R) 

CR“CMPLX( R,0 .0 ) 

ACMP«COS(GM*PHXP) 

CG2M» (CMPLX( ACMP,0 .0 ) *C EXP (CARG) -CMP LX (1.0, f.l) )/CR 

RETURN 

END 

SUBROUTINE  ELLPTC ( TL , TU,NGQ, CANS) 

IMPLICIT  COMPLEX  (C) 

REAL  COSG 

COMMON/TPI/TWOP I , MAX P , ER, P I 

COMMON 'FNC/RF, ZF ,TF, RSL ,ZSL, DEL, SING, COSG ,GM 

COMMON/SLF/ISELF 

EXTERNAL  CGME 

CALL  CGQl (CGME, TL,TU,NGQ, CANS!) 

CA,NS'CANS1*CMPLX(2 ,0/PI  ,9 .0  ) 

IF( ISELF  ,EQ.  0)  RETURN 

TLLN-8,0 

TULN-0.0 

IF(TF  .GT.  TL)  TLLN*  (TF-TL)  *ALOG (TF-TL) 

I P (TF  ,LT.  TU)  TULN= (TU-TF) ‘ALOG (TU-TF) 
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BR-13340 

BR- 1.3350 

BR-13360 

BR-I3370 

BR-13380 

BR-13390 

BR- 13400 

BR-13410 

BR-13420 

3R-13430 

BR-I3440 

BR-13450 

BR-13460 

BR-13470 

BR-13<90 

BR-13490 

BR- 13500 

BR-13510 

BR-13520 

BK-1J5J0 

Bk— 1 JS40 

BR-135S0 

B *13560 

BR-135/0 

BR-13580 

BR-13590 

Bk-UWK 

BR-13610 

BR-1362  0 

BR-1363  0 

BR-13M0 

BR-13650 

BR-13660 

BK-Ub/0 

BR-13680 

BR-Ub90 

BR-13700 

BR-U/10 

BR-13720 

BR-13730 

BR-13740 

BR-13/50 

BR-13760 

BR-13770 

3R-13780 

BR-13790 

BR-13800 

BU-13810 

BR-1JU20 

BR-13830 

UR-13840 

BR-13850 

BR-13860 

BR-13870 

BR-138  80 

BR-13890 

BR-13900 

BR-13910 

BR-13920 

BR-13930 


109 


mis  tm  is  jmi(#iMi^m<!W<asa 

UKjg  ^WY»V  niwiSHSD  TO  OPO 


ANS- ( (TO-TL) * (1.0 -A LOG (DEL) ) -TLLN-TULN) /RF/PI 
CANS-CANS+CMPLX ( ANS , 0 . 0 ) 

RETURN 

END 

COMPLEX  FUNCTION  CGME(TP) 

IMPLICIT  COMPLEX  (C) 

REAL  COSG 

COMMON/FNC/RF , ZF ,TF, RSL , ZSL ,DEL ,S ING , COSG ,GM 

COMMON/SLF/ISELF 

RS-RSL+TP* DEL* SING 

ZS»ZSL+TP« DEL* COSG 

ZFMZS2-  (ZF-ZS)  **2 

RFMRS2- (RF-RS)  **2 

RN-RFMRS2+ZFMZS2 

RD- (RF+RS) **2  +ZFMZS2 

R2-SQRT(RD) 

R3-SQRT(RN) 

AM1-RN/RD 
ANS-E  LI Cl K ( AMI ) 

ANS-ANS/R2 

IF(ISELF  ,EQ.  0)  GO  TO  10 
ANS-ANS+ALOG(R3) /2.0/RF 
10  CGME*CMPLX( ANS ,0 .0 ) 

RETURN 

END 

SUBROUTINE  ELLPTR(TL,TU,NGQ, CANS) 

IMPLICIT  COMPLEX  (C) 

REAL  COSG 

COMMON/TPI/TWOPI ,MAXP,ER,PI 

COMMON/FNC/RF , ZF , TF ,RSL , ZSL , DEL, SING , COSG ,GM 

COMMON/SLF/ISELF 

EXTERNAL  CGMER 

CALL  CGQ1 (CGMER, TL,TU,NGQ, CAN SI) 
CANS-CANS1*CMPLX (2 .0/PI ,0.0) 

IF(ISELF  . EQ.  0)  RETURN 

TLLN-0.0 

TULN-0.0 

IF(TF  .GT.  TL)  TLLN- (TF-TL)  *ALOG(TF-Ti.) 

IF(TF  .LT.  TU)  TULN«(TU-TF)*ALOG(XU-TP) 

ANS- ( (TU-TL) * (1.0-ALOG(DEL) ) -TLLN-TULN) /PI 
CANS-CANS+CMPLX( ANS,0.0) 

RETURN 

END 

COMPLEX  FUNCTION  CGMER(TP) 

IMPLICIT  COMPLEX  (C) 

REAL  COSG 

COMMON/FNC/R  F , Z F , TF , RSL , ZSL , D EL , S X NG  , COSG , GM 

COMMON/SLF/ISELF 

RS-RSL+TP* DEL* SING 

ZS-ZSL+TP* DEL* COSG 

ZFMZS2- (ZF-ZS) **2 

RFMRS2- (RF-RS) **2 

KN-RFMRS2+ZFMZS2 

PD-  (RF+RS)  “2+ZFMZS2 

R2-SQRT(RD) 

R3-SQRT(RN) 

AM1-RN/RD 

ANS-ELICIK(AMI) 

ANS-ANS*RS/R2 

IFdSELF  .EQ.  0)  GO  TO  10 

ANS»ANS+ALOG;R3)/2.0 
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BR-14260 
BR-14270 
BK-14  280 
BR-1429  0 
Bk-14300 
BR-14310 
BR- 14320 
BR-14330 
BR-14134# 
BR-14350 
BR-14360 
BR-14370 
BR-14360 
BR-14390 
BR-14400 
BR-14413 
BR-14420 
BR-14430 
BR-14440 
3R-14450 
BR- 14460 
BR-14470 
BK-14480 
BR-1449H 
BR-14500 
BR-14510 
BR-14520 
BR-14530 
BR-14540 
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10  CGMER-CMPLX (ANS,0.0 ) 

RETURN 

END 

COMPLEX  FUNCTION  CO POR( DUMMY) 

FUNCTION:  (1/(U0))*(D(G1+G2)/D(R0)) 

IMPLICIT  COMPLEX  (C) 

REAL  COSGS 

C0MM0N/FNC/RF,ZF,TF,R5L,ZSL, DEL, SINGS, CO SGS,GM 

COMMON/ FRQ/P I , AK 1 , AK2 , SL1 , SL2 , OMEGA 

COMMON/ARC/TP 

COMMON/ANG/PHIP 

RS=RSL+TP*DEL*S CNGS 

ZS=*ZSL+TP*DEL*COSGS 

R-SQRT (RF*RF+RS*RS+ (ZF-ZS) **2-2.0 *RF*RS*COS (PHIP) ) 

RK1=AK1*R 

RK2«AK2*R 

C ARGl “C  MP  LX ( 0 . 0 ,-RKl) 

CARG2*CMPLX(0.Ci  ,-RK2) 

C1PJK1-CMPLX  (1 ..  0 , RK1) 

C1PJK2«CMPLX(1.B,RK2) 

CB=C1?JK1*CEXP(CARG1) +C1PJK2*CEXP (CARG2) 

CGPOR— CB/CMPLX(R**3,0.0) 

RETURN 

END 

COMPLEX  FUNCTION  G0 (TP) 

IMPLICIT  COMPLEX  (C) 

C UMMON/  T Pl/TWOPI , MAX  P , E R,  P I 
COMMON/ A RC/T 
EXTERNAL  CGU 
T-TP 

CALL  TRPADP ( CG0 , 0 . 0 , P I , ER, MAXP , I ER, CANS ) 

G0-CANS/C  MPLX ( TWOPI ,0.0) 

RETURN 

END 

COMPLEX  FUNCTIUN  G1 (TP) 

IMPLICIT  COMPLEX  (C) 

COMMON/ I SWTCH/ISWR 

COMMON/TPI/TWOPI , MAXP , ER, PI 

COMMON/ARC/T 

EXTERNAL  CGI 

ISWR-0 

T-TP 

CALL  TRPADP (CGI, 0.0, PI, ER, MAXP , I ER, CANS) 

G1 -CAN S/C MPLX (TWO? 1,0 ,0) 

RETURN 

END 

COMPLEX  FUNCTION  GIR(TP) 

IMPLICIT  COMPLEX  <C) 

COMMON/iS«TCH/ISWR 

COMMON/TPI/TWOPI, MAXP, ER, PI 

COMMON/ARC/T 

EXTERNAL  CGI 

ISWR-1 

T-TP 

CALL  TRPADP (CGI, 0. 0,P I, ER, MAXP, IER, CANS) 

G'JR-CANS/C  MPLX  (TWOPI  ,0,0) 

RETURN 

END 

complex  function  G2  (TP) 
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IMPLICIT  COMPLEX  (C) 

COMMON/TPI/TWOPI ,MAXP,ER,PI 

COMMON/ARC/T 

COMMON/ I SWTCH/ISWR 

EXTERNAL  CG2 

T»TP 

ISWR*B 

CALL  TRPADP (CG2,0 ,0 ,PI,ER,MAXP,I EP ,CANS) 
G2»CANS/CMPLX(TWOPI,0.0) 

RETURN 

END 

COMPLEX  FUNCTION  C2 R ( TP) 

IMPLICIT  COMPLEX  (C) 

COMMON/TPI/TWOPI , MAXP , ER,  PI 

COMMON/ARC/T 

COMMON/ I SWTCH/ I SWR 

EXTERNAL  CG2 

T*TP 

ISWR*1 

^ ™PADP{CG2 , 0 . 0 , p 1 , ER, MAXP, IER,C ANS) 
G2R»CANS/CMPLX(TWOPI.0.0) 

RETURN 

END 

COMPLEX  FUNCTION  G3.'TP) 

IMPLICIT  COMPLEX  (C) 

COMMON/TPI/TWOPI , MAXP, ER, PI 

COMMON/ARC/T 

EXTERNAL  CG3 

T"TP 

Siicl»S™  :!!i  «*>»»».»••.<»» ) 

RETURN 

END 

COMPLEX  FUNCTION  G4R(TP) 

IMPLICIT  COMPLEX  (C) 

COMMON/TPI/TWOPI , MAXP, ER,PI 

COMMON/ARC/T 

COMMON/I SWTCR/ISWU 

EXTERNAL  CG4 

T*«TP 

ISWR-1 

RETURN 

END 

COMPLEX  FUNCTION  G5 (TP) 

IMPLICIT  COMPLEX  (C) 

COMMON/TPI/TWOPI, MAXP, ER,  PI 

COMMON/ARC/T 

C OMMON/ 1 SWTCH/ 1 SWR 

EXTERNAL  CG5 

T-TP 

ISwR»» 

1 Ca 5 • 0 • B ' p I » E R , MAX  P , I ER , C ANS ) 
G5«CANS/CMPLX(TWOPI,0,0)  ' 

RETURN 
END 

COMPLEX  FUNCTION  G5RJTP) 

IMPLICIT  COMPLEX  *C) 

COMMON/TPI/TWOPI, MAXP, ER, PI 
COMMON/ARC/T 
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COMMON/ ISWTCH/ 1 SNR 
EXTERNAL  CG5 
T-TP 
I8WR-1 

CALL  TRPADP(CG5 ,0 .0 ,PI ,ER, MAXP,IER,CANS) 
G5R-CANS/CMPLX(TWOPI,0 ,0 ) 

RETURN 

END 

COMPLEX  FUNCTION  G6{TP) 

IMPLICIT  COMPLEX  (C) 

COMMON/TPI/TWOPI , MAXP , ER, PI . 

COMMON/ARC/T 
EXTERNAL  CGb 
T«TP 

Call  trpadp ( cgs ,0 . 0 ,p 1 , er, maxp , 1 er,c ans ) 
Gb-CANS/CMPLX( TWOPI ,0 .0) 

RETURN 

END 

COMPLEX  FUNCTION  CUD (PHIP) 

IMPLICIT  COMPLEX  (C) 

REAL  COSGS 
COMMON/ANG/PHI 
COMMON/ARC/TP 
COMMON/SLF/ISELF 

COMMON/FNC/RF,ZF,TF,RSL,ZSL, DEL, SINGS, COSGS, GM 
COMMON/FRQ/PI , AX1 , AX  2 , SL1 , SL2 , OMEGA 
PHI-PHIP 

SMPSP-SIN(GM*PHIP) *SIN(PHIP) 
CINTG-CMPLX(SMPSP,0.0) *CGPOR(XXXXXX) 

IFtlSELF  .EQ.  0)  GO  TO  10 

ANUM-2.0*GM*PHIP*PHIP 

DEN- ( ( (TF-TP) *DEL) **2+(RF*PHIP)  **2) **1,5 

CINTG-CINTG+CMPLX (ANUM/DEN,0 .0 ) 

10  CGt-CINTG 
RETURN 
END 

COMPLEX  FUNCTION  CGl(PHJV) 

IMPLICIT  COMPLEX  (C) 

REAL  COSGS 
COMMON/ISWTCI  /ISWR 
COMMON/ANG/PHI 
COMMON/ARC/TP 

COMMON/FNC/RF,ZF,TF,RSL,ZSL, DEL, SINGS, COSG8,GM 
COMMON/FRQ/PI , AK1 , AK2 , SL1  , SL2  , OMBGA 
PHI-PHIP 

ZS-ZSL+TPwDEL*COSGS 
SMPSP-S1N(GM*PHIP) *SIN(PHIP) 

IFdSMR  .EQ,  1)  SMPSP-SMP5P"  (RSL+TP*DEL*SINGS) 
CG1-CMPLX( (ZF-ZS) *SMPSP,0 .0 ) »CGPOR(XXXXXX) 
RETURN 
END 

COMPLEX  FUNCTION  CG2(PHIP) 

IMPLICIT  COMPLEX  (C) 

REAL  COSGS, CMPCP 
COMMON/ANG/PHI 
COMMON/ ISWTCH/ ISWR 
COMMON/ARC/TP 
COMMON/SLF/ISELF 

COMMON/FNC/RF, ZF,TP,RSL , ZSL , DEL, SINGS, COSGS ,GM 
COMMON/ F RQ/ P I , AK1 , AK  2 , SL1 , SL2 ,OHBGA 
PHI-PHIP 


BR-15//0 

BR -lb /HR 

BK-ib/90 

BK-J.5800 

BR-15810 

Bft-15820 

BK-lbUJH 

NM-13040 

BK-lbbb* 

BR-15860 

BR-15870 

BR-15000 

BR-15890 

BK-.IS900 

BR-15810 

BR-15920 

BR-15930 

BK-13940 

BR-15950 

BR-159b0 

BM-idy/i 

8R-15980 

BR-15990 

BR-16000 

BR- 16010 

BR-16020 

BR-16030 

BR-16040 

BR-16050 

BR-16060 

BR-lbV/0 

BR-16080 

BR-16090 

BR-16100 

BR-lbl 10 

BR-lbl 20 

BR-161 30 

BR-16140 

BR-lbl 50 

BR-16160 

BR-16170 

BR-lblH  0 

BR-1619M 

BR-16200 

BR-16210 

BR-16220 

BR-16230 

BR-16240 

BR-16250 

BR-lb2t0 

BR-lb2/0 

BR-16280 

BR-1629  0 

BR-16300 

BR-16310 

BR-16320 

BR-16330 

BR-16340 

BP- 16350 

BR-16360 

BR-16370 
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ZS-ZSL+TP*DEL*COSGS 
CMPCP-COS(GM*PHIP) *CCS(PBIP) 

IF(ISWR  .EQ.  1)  CMPCP-CMPCP* (RSL+TP*DEL*8IMGS) 
CINTG-CMPLX(  (ZP-ZS)  *CMPCP,0.0)  *CGPOR(X XXXXX) 
IFflSELF  .EQ.  0)  GO  TO  II 
ANOM-2.0* (TF-TP) *DEL*COSGS 
IP(ISWR  .EQ.  1)  ANUM-ANUM*RF 
DEN-  ( { (TF-TP)  *DEL)  **2-*(RF*PHIP)  **2)  **1,5 
. CINTG-CINTG+CMPIK  (ANUM/DEN,0 .0 ) 

10  CG2-CINTG 
RETURN 
END 

COMPLEX  FUNCTION  CG3(PHIP) 

IMPLICIT  COMPLEX  (C) 

REAL  COSGS 
COMMON/ANG/PRI 
COMMON/ARC/TP 
COMMON/SLF/ISELF 

COMMON/FNC/RF, ZF,TF,RSL, ZSL, DEL, SINGS, COSGS,GH 
COMMON/FRQ/PI ,AK1 ,AK2,SL1 ,SL2 , OMEGA 
PHI-PHIP 

RS-RSL+TP* DEL* SINGS 
SMPSP-SIN(GH*PHIP)  *SIN(P!!IP) 

C INTG-C MP LX ( RS* S MP SP , 0 . 0 ) *CG?OR (X XXXXX) 

IF(ISELF  .EQ.  0)  GO  TO  10 

ANUM-2.0*RP*GM*PHIP*PHIP 

DEN- ( ( (TF-TP) *DEL) **2+(RF*PHIP) **2) **1.5 

C INTG-C INTG+CMP LX (ANUM/DEN, 0.0) 

10  CG3-CINTG 
RETURN 
END 

COMPLEX  FUNCTION  CG4(PHIP) 

IMPLICIT  COMPLEX  (C) 

REAL  COSGS, CMP 
COMMON/ ANG/PH I 
COMMON/ I SWTCH/ I oWR 
COMKON/ARC/TP 
COMMON/SLF/ISELF 

COMMON/FNC/RF , Z F ,TF ,RSL , ZSL , DEL, SINGS , COSGS, GM 
COMMON/FRQ/P I , AK1 , AK  2 , SLl , SL2 , OMEGA 
PHI-PHIP 

RS-RSL+TP"DEL»SI«GS 
CMP-COS (GM*PHIP) 

IFdSWR  .EQ.  1)  CMP-CMP *RS 

C INTG-CMPLX ( ( RF-RS) *CMP , 0 . 0 ) *CGPOR (X  XXXXX) 

IFdSELP  .EQ.  0)  GO  TO  10 

ANUM-2.0* (TF-TP) *DEL*SINGS 

IFdSWR  .EQ.  1)  ANUM-ANUM*RF 

DEN- ( ( (TF-TP) *DEL) **2+(RF*PHIP) **2)  **1.5 

C INTG-C INTG+CMP LX (AHUM/DEN, 0.0 ) 

10  CG4-CINTG 
RETURN 
END 

COMPLEX  FUNCTION  CG5(PHIP) 

IMPLICIT  COMPLEX  (C) 

REAL  COSGS 
COMMON/ I SWTCH/ISWR 
COMMON/ ANG/PH I 
COMMON/ARC/TP 
COMMON/SLF/ISELF 

COMMON/FNC/RF, ZF,TF,RSL, ZSL, DEL, SINGS, COSGS ,GM 


BR-16380 
BR-16390 
BR- 16400 
BR-16410 
BR-16420 
BR-16430 
BR-16440 
BR-16450 
BR- 16460 
BR-16470 
BR-16480 
BR-16490 
BR-16500 
BR-16510 
BR-16520 
BR-16530 
BR-16540 
BR-16550 
BR-16560 
BR-16570 
BK-16500 
BR-16590 
BR-16600 
BR-16610 
BR-16621 
BR- 16630 
BR-16640 
BR-16650 
BR-16660 
BR- 16670 
BR-16600 
BR-16690 
BR-16700 
BR-16/10 
BR- 16/20 
BR-16/30 
BR-16/40 
BR-16/50 
BR-14/60 
BR-16770 
BR-167B0 
ER-16790 
BR-16000 
BR*  16010 
BR-16020 
BK- 16830 
BR-16840 
BR*  16850 
BR-16869 
PR-16870 
8R-16MN0 
BR-16890 
BR-.'  690  0 
BR-} 6910 
BR-.6920 
BR-16930 
ER-1694  0 
BR-16550 
Bit- 16960 
BR-16570 
BR- 16980 
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COMMON/FRQ/PI,AKl,AK2,SLl,SL2,OMBGA 

PHI-PHIP 

S2PCM?=COS (GM*PHIP) * (S IN (PHIP/2 .0 ) **?) 

IF( ISWR  .EQ.  1)  S2PCMP-S2PCMP* (RSL+TP*DEL*SINGS) 

C INTG-CMPLX ( S2  PC  MP , 0 . 0 ) *CGPOR ( XXXXXX ) 

IF(ISELF  ,EQ.  0)  GO  TO  10 

ANUM-PHIP*PHIP/2 .0 

IF (ISWR  .EQ.  1)  ANUM-ANUM*RF 

DEN- ( ( (TF-TP) *DEL) **2+(RF*PHIP) **2) **1.5 

CINTG-C INTG+CMPUt  (ANUM/DEN , 0.0 ) 

10  CG5-CINTG 
RETURN 
END 

COMPLEX  FUNCTION  CG6(PHIP) 

IMPLICIT  COMPLEX  (C) 

REAL  COSGS ,CMFCP 
COMMON/ANG/PHI 
COMMON/ARC/TP 
COMMON/S LF/ I SELF 

COMMON,/ FNC/RF , ZF,TF ,RSL , ZSL , DbL,  SINGS,  COSGS ,GM 
COMMON/ FRQ/ PI , AK1, AK2,SL1 , SL2 , OMEGA 
PHI-PHIP 

RS-RSL+TP* DEL* SINGS 

cmpcp-cos;gm*phip) *cos(phip) 

C1NTG»CMPLX( (RF-RS) *CMPCF ,0 .0 ) *CGPOR(XXXXXX) 

IF(ISELF  .EQ.  0)  GO  TO  10 

ANUM-2.0* (TF-TP) *DEL*SINGS 

DEN- ( ( (TF-TP) *DEL) **2+( RF*PHIP)  **2)  **1.5 

CINTG-CINTG+CMPLX ( ANUM/DEN ,0 .0 ) 

10  CGb-CINTG 
RETURN 
END 

COMPLEX  FUNCTION  CANALY(TL,TU,I ,ISWR) 

ROUTINE  TO  ADD  RESULT  OF  ANALYTICAL  INTEGRATION  OF  THE 

SINGULAR  PORTION  OF  THE  INTEGRANDS 

COMMON/FNC/RP,ZF,TF,RSL, ZSL, DEL, SINGS, ACOSGS.GM 

COMMON/FRQ/PI ,AK1 , AK2,SL>1 ,SL2,OMBGA 

ID-I+1 

TUMTDT- (TU-TF) *DEL 
TLMTDT- (TL-TF) *DEL 
RPI-RF*PI 

IF(TUMTDf)  10,  10,  20 
10  TUI-0.0 
TU2-0.0 
GO  TO  30 

20  RAD-SQRT (TUMTDT*  *2+RPI*  *2 ) 

TUl-ALOG ( RPI+RAD) 

TU2-ALOG (TUMTDT) 

30  IF (TLMTDT)  50,  40,  40 
40  TL1-0.0 
TL2-0.0 
GO  TO  60 

50  RAD-SQRT (TLMTDT* *2  +RPI * *2 ) 

TLl=ALOG(  RADi-RPI) 

TL2-ALOG! -TLMTDT) 

r<0  TERM-  (TU-'.'F)  *(TUl-TU2)  + ( TF-TL)  * (TL1-TL2 ) 

TERM-  -TERM/ (PI*RF) 

IF ( ISWR  .EQ.  0)  TERM-TERM/RF 

GO  TO  (70,  100,  100,  80,  100,  90)  ID 


BR-16990 
BR-17000 
BR-17010 
BR-17020 
BR-17030 
BR-1 7040 
BR-I7050 
BR-17060 
BR-I7070 
BR-17080 
BR-1/ 090 
BR-17100 
BR-1 7110 
BR-17120 
BR-17130 
BR-17140 
BR-17150 
BR-17160 
BR-17170 
BR-1/18  0 
BR-1/ 19  0 
BK-1/200 
BR-1 ,210 
BP.-17  220 
BR-17230 
BR-17240 
BR-17250 
BR-17260 
BR-1727  0 
BR-1 7 28  0 
BR-1/ 29  0 
BR-1/ 300 
BR-17310 
BR-17320 
BR-17330 
BR-17340 
BR-17350 
BR-1 73 60 
BR-1737  0 
BR-1 7 380 
BR-17390 
BR-17400 
BR-1  J4  111 
BR- 17420 
BR-17430 
BR-17440 
BR-17450 
BR-17460 
BR-17470 
BR-17480 
BR-17490 
BR-1/ 500 
BR-17510 
BR-1/520 
BR-17530 
BR-1 7 546 
BR-17550 
BR-175S0 
BR-17570 
BR-17580 
BR-17590 
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SUBROUTINE  CSMINV ( A ,NDIM ,N ,DETERM ,COND ,IERR)  HI- 

****************************  » ***************************************  ♦•f'm- 

CSMINV  IS  A SUBROUTINE  WHICH  WILL  ACCEPT  A SINGLE  PRECISION  MI- 
COMPLEX  MATRIX  AND  RETURN  THE  INVERSE  OF  THE  MATRIX  IN  ITS  MI- 
PLACE.  THE  SUBROUTINE  WILL  ALSO  COMPUTE  THE  NORMALISED  ■ M> 
DETERMINANT  OF  THE  MATRIX,  AND  THE  INVERSE  CONDITION  NUHBER  MI- 
OF  THE  MATRIX.  MI- 

MI- 

N - THE  ORDER  OF  THE  MATRIX  TO  BE  INVERTED  MI- 

A - COMPLEX  DOUBLE  PRECISION  INPUT  MATRIX  (DESTROYED)  MI- 
THE  INVERSE  OF  A IS  RETURNED  IN  ITS  PLACE.  MX- 

NDIM  - THE  SIZE  TO  WHICH  A IS  DIMENSIONED  IN  THE  CALLING  MI- 
PROGRAM  mi- 

determ  - THE  NORMALIZED  DETERMINANT  OF  A WHICH  IS  RETURNEDMI- 
COND'-  THE  INVERSE  OF  MITTRA'S  CONDITION  NUMBER  OF  MI- 

THE  MATRIX.  MI- 

IERR  - ERROR  INDICATOR  WHOSE  VALUE  IS  ZERO  UNLESS  TOO  MI- 
LARGE  A MATRIX  IS  PASSED  ( .GT,  258  X 258),  MI- 

IN  WHICH  CASE  IERR-1  MX- 

MI- 

PREPARED  BY  MICHAEL  G.  HARRISON  E.E.  DEPT  JUNE  23,  1972  MI- 

MI- 

********* ********* ************** ****** ************* *•••****•**••• •••••■nx- 

complex  A(NDIH,NDIM) ,PIVOT(250)  ,AMAX,T, SWAP, DETERM, U MI- 

INT£GER*4  IPIVOT(250) ,INDEX(250,2i  MI- 

REAL  TEMP, ALPHA (2 58)  MI- 

MI- 

INITIAL1ZATION  MI- 

MI- 

IERS»0  MI- 

IF{NDIM.LE.250)  GO  TO  5 NI- 

IERR-1  MI- 

WRITE{5,4)  NDIM  MI- 

4 FORMAT ( '0CSMINV  ERROR.  ATTEMPT  TO  INVERT  A MATRIX  '14,  MI- 

1*  ON  A SIDE,’/'  WHEN  250  X 250  IS  THE  MAXIMUM  ALLOWED.')  MI- 

RETURN  MI- 

5 CONTINUE  M i- 

DETERM  « CMPLX (1.0, 0.0)  MI- 

SUHAXA-0.  Mi- 

DO  20  J=.l,N  Ml- 

ALPHA (J ) *0.0  MI- 

SUMROW»0.  MI- 

DO  10  I»1,N  MI- 

ALPHA(J) =ALPHA(J) +A(J,I) * CONJG ( A (J , I) ) MI- 

10  SUMROW=SUMROW  + CABS  ( A (J,  I))  MI- 

ALPHA(J) » SQRT(ALPHA(J) ) Mi- 


ll 

20 

30 

48 

50 

61 

71 

80 

31 
108 
119 
128 
138 
148 
158 
168 
178 
188 
198 
288 
218 
228 
238 
248 
258 
268 
278 
280 
298 
388 
318 
328 
338 
340 
358 
368 
378 
388 
398 
400 
410 
428 
430 
4 40 
450 
460 
470 


116 


ror.  o o o n n r>  n non  n r.  n 


tB18  MOB  IS  BIST  QUALITY  PRACTICABLE 
iROi  oopy  juwishto  to  roc 


60 

80 

85 

180 

185 

140 

200 

260 


3 30 

350 

?80 

400 

453 

550 


T (SUMROK.GT.SUMAXA)  SUMAXA-SUMROW 

MI- 

480 

1PIV0T(J) *0 

KI- 

4 90 

DO  600  1*1, N 

MI- 

500 

Ml- 

610 

SEARCH  FOR  PIVOT  ELEMENT 

MI- 

520 

Ml- 

610 

AMAX=CMPLX(0. 0,0.0) 

MI- 

540 

DO  105  J-l  ,H 

MI- 

550 

IF  (IPIVOT(J) -1)  60,  105,  60 

MI- 

560 

DO  100  K»1,K 

MI- 

570 

IF  (IPIVOT(K) -1)  80,  100,  740 

MI- 

580 

TEMP-AMAX*  C ON JG  < AMAX ) -A ( J , K ) * CON JG ( A ( J , K ) ) 

MI- 

590 

IF (TEMP) 85,85,100 

MI- 

600 

IROW-J 

MI- 

610 

ICOLUM-K 

MI- 

620 

AMAX-A(J,K) 

MI- 

630 

CONTINUE 

MI- 

640 

CONTINUE 

MI- 

650 

I P I VOT ( I COLUM ) -I P IVOT ( I COLUM ) f 1 

MI- 

660 

MI- 

670 

INTERCHANGE  ROWS  TO  PUT  PIVOT  ELEMENT  ON  DIAGONAL 

MI- 

680 

MI- 

690 

IF  (I ROW- I COLUM)  140,  260,  140 

MI- 

700 

DETERM— DETERM 

MI- 

710 

DO  200  L»1,N 

HI- 

720 

SWAP* A ( I ROW, L) 

MI- 

730 

A ( IROW,L) “A ( I COLUM, L) 

MI- 

740 

A ( ICOLUM ,L) -SWAP 

MI- 

750 

SWAP-ALPHA ( IROW) 

MI- 

760 

ALPHA (I ROW) “ALPHA{ ICOLUM) 

MI- 

770 

ALPHA(.tCOLMM)  -SWAP 

MI- 

780 

INDEX  (1,1)  -IF.OW 

MI- 

790 

INDEX (I, 2) -ICOLUM 

HI- 

800 

PIVOT ( I ) -A (ICOLUM , ICOLUM) 

MI- 

810 

U - PIVOT(I) 

MI- 

820 

ALPHAI- ALCHA (ICOLUM) 

MI- 

830 

MI- 

840 

• THE  FOLLOWING  SUBROUTINE  CALL  IS  POR  UNDERFLOW  PROTECTION 

DURING  MI- 

850 

CALCULATION  OF  THE  NORMALIZED  DETERMINANT 

MI- 

860 

MI- 

870 

CALL  DTRMNT(DETERM,U,ALPHAI) 

MI- 

880 

TEMP-PIVOT(I)  * CONJG(PIVOTd)  ) 

MI- 

890 

IF (TEMP) 330,720,330 

MI- 

900 

MI- 

910 

DIVIDE  PIVOT  ROW  BY  PIVOT  ELEMENT 

MI- 

920 

MI- 

930 

A (ICOLUM,  ICOLUM)  - CMPLX  (1 .0  ,0 .0  )' 

MI- 

940 

DO  350  L»1,N 

MI- 

950 

U - PIVOT(I) 

MI- 

960 

A (ICOLUM ,L)  - A (ICOLUM, L)/U 

MI- 

970 

MI- 

980 

REDUCE  NON- PIVOT  ROWS 

Ml- 

990 

MI- 

1000 

DO  558  Ll-1 ,N 

MI- 

1010 

IF(Ll-ICOLUM)  400,  558,  400 

MI- 

1020 

T-A(Ll, ICOLUM) 

MI- 

1030 

A(L1, ICOLUM) - CMPLX(0. 0,0,0) 

MI- 

1040 

DO  450  L- 1 ,N 

Ml- 

1860 

V - A ( ICOLUM, L) 

MI- 

1060 

A (LI , L)  - A(L1,L)-U*T 

MI~ 

10/0 

CONTINUE 

MI- 

1080 
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600  CONTINUE  MI- 

MI- 

INTERCHANGE  COLUMNS  MI- 

MI- 

620  DO  710  1*1, N MI- 

L-N+l-I  MI- 

IP  (INDEX(L,1)-INDEX(L,2) ) 630,  710,  630  MI- 

630  JROW*INDEX(L,l)  MI- 

JCOLUM* INDEX (L, 2}  MI- 

DO  705  K*1,N  MI- 

SWAP*  A (K,J  ROW)  MI- 

A(K,JROW)  *A(K,JCOLUM)  MI- 

A (K , J COLUM) *SWAP  MI- 

705  CONTINUE  MI- 

710  CONTINUE  MI- 

SUMAXI*0.  MI- 

DO  910  1-1  ,N  MI- 

SUMROW-0.  MI- 

DO  900  J-l ,N  MI- 

900  SUMROW=SUMROW  + CABS(A(I,J))  MI- 

IF(SUMROW.GT.SUMAXI)  SUMAXI-SUMROW  MI- 

910  CONTINUE  MI- 

COND  - 1./ (SUMAXA*SUMAXi ) MI- 

RETURN  MI- 

720  WRITE(5 ,7  30)  MI- 

730  FORMAT  ( ’0  * ,10(  •*******■•)/  '0  MATRIX  IS  SINGULAR’ /'0*  ,10  (•*******') ) MI- 
740  RETURN  MI- 

END  MI- 

SUBROUTINE  DTRMNT(DETERM, U, A)  MI- 

£**********•************************** **»*••*#**•••»*«*•*•»••*••*«*• •••*)(  J_ 

C HI- 

C THE  SOLE  PURPOSE  OF  THIS  ROUTINE  IS  TO  SCALE  THE  NORMALIZED  KI- 

C DETERMINANT  SO  THAT  MACHINE  UNDERFLOWS  WILL  NOT  OCCUR.  MI- 

C THIS  IS  NECESSARY  BECAUSE  THE  DYNAMIC  RANGE  OF  THE  DEC-1177  MI- 

C MACHINE  IS  ONLY  ABOUT  10* *-38  TO  10**+39.  MI- 

C MI- 

C THE  PARAMETER  ISCALE  MUST  BE  RETURNED  VIA  COMMON  TO  ANY  MI- 

C PROGRAM  NEEDING  THE  VALUE  OF  THE  NORMALIZED  DETERMINANT.  MI- 

C MI- 

C THE  VALUE  OF  THE  NORMALIZED  DETERMINANT  IS  THE  VALUE  RETURNED  MI- 

C BY  CSMINV  TIMES  TEN  RAISED  TO  THE  POWER  (-1 SCALE*TEN) . MI- 

C Mi- 
ll MI- 

C IF  CSMINV  IS  CALLED  MORS  THAN  ONOE  BY  A PROGRAM,  "ISCALE"  MUST  BE  MI- 

C INITIALIZED  TO  ZERO  FOR  EACH  CALL  AFTER  THE  FIRST.  MI- 

C MI- 

(jxt**********************************  ************** *********************  mi- 

complex  DETERM, U MI- 

COMMON/SCAFAC/I SCALE  MI- 

DATA  I SCALE/0/  MI- 

IF (CABS (DETERM)  0GT.  l.E-10)  GO  TO  100  MI- 

DETERM*DETERM*1 . E10  MI- 

ISCALE* I SCALE+1  MI- 

100  DETERM=DETERM*U/CMPLX( A ,0 ,0 ) MI- 

RETURN  MI- 

END  MI- 


1090 

1100 

1110 

1120 

1130 

1140 

1150 

1160 

1170 

1180 

1190 

1200 

1210 

1220 

1230 

1240 

1250 

1260 

1270 

1280 

1290 

1300 

1310 

1320 

1330 

1340 

1350 

1360 

1370 

1380 

1390 

1400 

1410 

1420 

1430 

1440 

1450 

1460 

1470 

1480 

1490 

1500 

1519 

1520 
1530 
1540 
1550 
1560 
1570 
1580 
1590 
1600 
1610 
1620 
16  38 
1640 
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SUBROUTINE  CGQ1 (CF,XL, XU,N, CVAL) 

«**,!********!>***************************,►******»**»*****»*****»***■»** 

PERFORMS  INTEGRATION  OF  A FUNCTION  OF  A SINGLE  VARIABLE 
BY  GAUSSIAN  QUADRATURE. 

N - ORDER  OF  GUASSIAN  QUADRATURE  APPROXIMATION 
(2,4,8,10/12,16,32) 

CF  =*  EXTERNALLY  SUPPLIED  FUNCTION ...  .MUST  BE  FUNCTION  OF  ONE 
VARIABLE  FOR  CQG1 . 

XL  ■»  LOWER  BOUND  OF  VARIABLE 
XU  = UPPER  BOUND  OF  VARIABLE 
CVAL  - RESULTING  VALUE  OF  THE  INTEGRATION 


CF  MUST  BE  LISTED  IN  AN  EXTERNAL  STATEMENT 
PREPARED  BY  MICHAEL  G.  HARRISON  E.E.  DEPT 


JUNE  22,  1972 


Q 

**Q 

Q 

Q 

Q 

Q 

Q 

Q 

Q 

Q 

Q 

Q 

Q 

Q 

Q 

Q 

Q 

Q 

Q 


IMPLICIT  COMPLEX  (C)  Q 

DIMENSION  Ql (52) ,Q2(24) ,Q3<32) ,NQ(8) ,NS(8) ,QG{108)  Q 

EQUIVALENCE  (Ql( 1 ) ,QG(1 ) ) . (02(1 ) ,QG( 5 3) ) , (Q3 (1 ) ,QG( 77) ) Q 

DATA  Ql  Q 

$ /.283675134594812B82E0,0.5E0, .4 3056815 57 9702629E0 , Q 

5. 17392742256872693E0, .1 699905217924281 3E0 , . 3 2607257743127307E0,  Q 
? 0.4 801445 28 248 768 12E0 , .5 061 426 8 14518 81 30E-1, .3 983332 3870 68 1337E0,Q 
$.1 1119051722668724E0,. 26276620495 816449E0,. 156853322938 94364E0,  Q 
$. 91717321 24782490 E-l,.  1813418  91689180 99E0,  .4 8695  3264  25  858586E0 , Q 

$. 3333567215434407E-1,. 43253 168334 4492 25E0,. 7 4725674 5752903E-1,  Q 
?. 3 39704784149  6122E0  ,.10954 3181 25 799102E0,  . 216697 69  70  646236E0,  Q 
5. 1 3463335915499818E0, . 7443716 9490 8l5605E-l,.l 4776211235737 64 4E0,  Q 
$0, 4 90780317 1 2335963 EC,. 2 358 76 68 19 3 25 591 4 E-l, .452058  6281852374  3E0 , Q 
$. 534696629976592 15E-1,. 384951337097152 34E0, .8 00391642716731 IE- 1,  Q 
S.29365897714330872E0 ,.l 0158371 336 153 296E0, .18391 57494990901 0E0,  Q 
S. 11674626826917 740E0, .6 261670 42557 34458E-1, .12457352 2906701 39E0,  Q 
$.4  9<I700467  49582  497E0  . .1357622970 58 770 47E-1,  .4 7228  751 1536616 29 E0,  Q 
$. 31126761969323946E-1,. 4 3281 560119391 587E0, .4 7579255841246392E-1,  Q 
$. 37 7 702 20 417750152 E0,. 6231448 56 27 76 69 36 E-l, . 3 08338 1 2220 13 21 8 7E0,  Q 
$,7479799440828837E~1 , ,229008  3888  2861 369E0 , ,8  45  78259697  5 0127 E-l, 

$. 14080177538 96 2946 E0, .9 13017075 2 246 17 9E-1, .47506254918818720E-1, 
$.947253052275 34 25 E-l  / 

DATA  Q2 

$ / 0.497 59 3609 998510 68 E+0 

* 0.4873642779856 54 75 E+0 

* 0.46913727 60 0136638 E+0 

* 0.443207 76350220 052 E+0 

* 0,4 10000992986951 46 E+0 

* 0,370062095 78 9 27 7 18 E+0 

* 0 .3  240468  2596  84 8778 E+0 

* 0.272710 /356 94 41 977E+0 

* 0,2168 96 75 38 1302 25 /E+0 

* 0.1 57521339 04 8081 69 5+0 

* 0, 955594 337368081 50E-1 

* 0.32028446431302813E-1 


,Q 

,Q 


DATA  Q3 


/0. 4 986319309247407  8E+0 
0,4 9 280 57 5 57 7 263  417  E+0  , 
0.48238112779375322E+0  , 

0.4 67 4 53 037 96 8 86 9 84 E+0  , 
8,448160577883026062+0  , 


Q 
Q 
Q 
Q 

,0.6 1706148 999935998E-2  ,Q 
0. 1 4265694  31 4 4668  3 2E- 1 ,Q 
0.22138719408 70990 3E-1 
0. 2 96 4929 2 457 7 18 89 0E-1 
0. 3 66 7 3 24 07 055401 5 3E-1  ,Q 
0.4  30950807 659 76638E- 1 ,Q 
0.4 8 3093 26 0520 5694 4E-1  ,Q 
0.5 37 2 21 35 05 79 8 281 7 E-l  ,Q 
0.5  775283  4026 8628 01 E-l  ,Q 
0.608352 36463901696E -1  ,Q 
0.6  29 18 7281 73 414 14 82- 1 ,Q 
0.6 396 90 97 67 3 37 60 7BE-1  /Q 

Q 

0.3  5093 050 04735 0483 E- 2 ,Q 
018 13719  736  545283  50E-2  ,Q 
0,1 26 960 326 546 3 10 30 E-l  ,Q 
0.17136931456510717E-1  ,Q 
0 .2  141794901 .1113  340E-1  ,Q 


10 

20 

30 

40 

50 

60 

70 

80 

90 

100 

110 

120 

130 

140 

150 

160 

170 

180 

190 

200 

210 

220 

230 

240 

250 

260 

270 

280 

290 

300 

310 

320 

330 

340 

350 

360 

370 

380 

390 

400 

410 

4 20 
431 
440 
450 
460 
470 
480 
490 
500 
510 

5 20 
530 
540 
550 
560 
570 
580 
590 
600 
610 


119 


orjnnnnrrnn 


*ii8  nn 
mo a 


* 0. 4 246 838 06 86 6 284 9 9E+0 

* 0 3 972 418 979 8 3971 20E+0 

* 0 .3 66091 05937 01 4484 E+0 

* 0 .33152  2133465107  60E+0 

* 0.29385787862038116  E+0 

* 0.253449 95 4 46 61 14 70 E+0 

* 0 .210675638065317  67E+0 

* 0 .16 59343011410638 2E+0 

* 0.1 1964 36 81 12606 8 5 4E+0 

* 0.7 223 59 807 9139 825 0E-1 

* 0 .2 415383284 38 69158E-1 
. $,NQ/2,4,8,10,12,16,24,32/, 

$ NS/1,3,7,15,25,37,53,77/ 

DO  300  L-1,8 

IF(N,EQ.NQ(L) ) GO  TO  301 

300  CONTINUE 
WRITE(5,905)  N 

905  FORMAT ( ' 0 CALLING  PARAMETER 
RETURN 

301  CONTINUE 
NP-NS(L) 

NE-NP+N-1 
AX-0.5* (XU+XL) 

BX-XU-XL 

CVAL-(0.,0.) 

DO  350  J-NP,NE,2 
DX-QG(J) *BX 

CVAL-CVAL+QGtJ+1)  *(CF(AX+DX)  +CF{AX-I 
350  CONTINUE 

CVAL-CVAL*BX 

RETURN 

END 


0.2 54 990 29 631 18 80 88E-1  ,Q  - 620 
0,293420 46 7392 6777 4 E-l  ,Q  - 630 
0.3 291 1111 3 8 81 8092 3E-1  ,Q  - 640 

0 .361728970  54 42 4253 E-l  ,Q  - 650 

0.3 909694789 353 5153E-1  ,Q  - 660 

0.4 1655 96 2 113 47 3 37 8E-1  ,Q  - 670 

0.4 38260 46 502 201906E-1  ,Q  - 680 

0.45586939347 8819  42E-1  ,Q  - 690 

0.46922199540402283E-1  ,Q  - 700 

0.47819360039637430E-1  ,Q  - 710 

0. 4 8 27004425 736390 0E-1  /Q  - 720 

Q - 730 

Q - 740 

U - /SB 
Q - 760 

U - / /0 

Q - 780 

INTEGRATION  NOT  POSSIBLE'//)Q  - 790 

Q - 808 

Q - 810 

U ~ o i* 
U - oJ0 
U - 840 

Q - 85  0 

U - oob 
Q - 870 

0 “ BOB 

) Q - 890 

Q - 9BB 
Q - 910 

Q - 92B 

U - 8 JB 


SUBROUTINE  CGQ1T(CF,XL,XU,N,CVAL,CVALR) 


SUBROUTINE  CGQ1T  IS  A VERSION  OF  THE  GAUSSIAN  QUADRATURE 
INTEGRATION  ROUTINE  "CGQ1"  WHICH  HAS  BEEN  MODIFIED  TO 
SIMULTANEOUSLY  INTEGRATE  THE  FUNCTIONS,  (CF)  AND  ( (RHO-PRIME) *CF) . 

THE  ARGUMENTS  ARE  THE  SAME  EXCEPT  FOR  THE  INCLUSION  OF  THE 
ADDITIONAL  RESULT  "CVALR." 


IMPLICIT  COMPLEX  (C) 

DIMENSION  Ql(52) ,Q2(24)  ,Q3 (3  2) ,NQ(8) ,NS{8) ,QG(108) 
REAL  COSGS 


QR- 

*QR- 

QR- 

U»- 

QR- 

QR- 

QR- 

QR- 

QR- 

QR- 

*QR- 

QR~ 

QR- 

QR- 


COMMON/FNC/RF, ZF,TF ,RSL , ZSL, DEL, SINGS , COSGS,GM  QR- 

EQUIVALENCE  (Q 1 (1 ) ,QG( 1 ) ) , (Q 2 (1 ) ,QG( 5 3) ) , (Q3 (1 ) ,QG{ 77) ) QR- 

DATA  Q1  QR- 

$ /, 2 88675 1345 94 81 28 8 2E0 ,0.5E0, ,4  30 5681 5 57 97 02 62 9 E0,  QR- 

$.1 7392742256872693E0,. 16999052 17924281 3E0,. 3 26072577431 27307E0,  QR- 

$ 0.4 80144928248 76812E0  ,. 506142681 4518 8130E-1,. 3 9833323870 681 337E0,QR- 


$.11 1190517226687 24E0 ,,26276620  495816449E0,  .15685332293394 364 E0 , QR- 
$. 9171732124782 490E-1,. 181341 89168918 099E0 , .4 869532642585858 6E0,  QR- 
$. 33335672154344 07 E-l,. 432531683344492 25E0  , .747256745752903E-1,  QR- 


$.3 397047841 'i9 61 22 E0,.l 0954 31 81 25 79 91 02 E0  , .2166976P70646236E0,  QR- 
$. 1 3463335915499818 E0, . 74437169490 815605E-I,. 14776 >11235737 64  4E0,  QR- 


10 
20 
3B 
4B 
50 
60 
70 
80 
90 
100 
110 
120 
13B 
140 
150 
160 
170 
180 
190 
200 
210 
2 20 
230 
240 
250 


120 


ttxs  nm  u as*  Qtuiijy  mmsASii 

f*tm  OOW-IWWtp^Tft^DC 


ft 


300 
905 

301 


$0. 49078031712335963E0,. 23587668193 25 5914E 

-1,  ,4 5 2058628 18 52 37 43E0, 

QR- 

260 

$. 5 3 46966 29976592 15 E- 1 , .3 84951 33709715234E0  , .8003916427 167311E-1, 

QR- 

270 

$.2936 58 97714330 872E0  10158 371336 15 3296 E0 

, .1839157494990 90 10 E0, 

QR- 

280 

$.116746 268 26917 74 0E0 , .6 26 1670 425 57 3 445 8E- 

1, .12457352 290670139E0, 

QR- 

290 

$. 4947004674958249 7E0 , .1 3576229705877047E- 

1,  .47228  7511!: 3661629E0, 

QR- 

300 

$. 3 1126 7619693 2394 6E-1, .4 3281 56011 9391 587E0 , .47579255041246392E-1. 

QR- 

310 

$.37770220417750152E0,.62314485627766936E- 

1 , .30893  81 222013  218 7E0, 

QR- 

320 

$.7 4 7979 9440 828 83 7E-1, . 2 290 08 3 8 88 28 61 36 9 E0 

, . 8457825969 750127E-1, 

QR- 

330 

$. 1408017753896 294 6E0, .9130170752246179E-1 

,. 47506254918818 72 0E-1, 

QR- 

340 

$.9 472530 522753 425E-1  / 

UR- 

350 

DATA  Q2 

UK- 

3 60 

$ 

/ 0.4 9 7593 6 0999 8 510 68 E+0 

, 0.6 17061 48 999935998E-2 

,QR- 

3"»0 

0. 4 37 36 427 7985 6b 4/ 5E+0  , 

0. 1426569431446 6832B-1 

»QR- 

380 

0.46913727600136638E+0  , 

0.2 2138 7194 08709903E-1 

,QR- 

390 

0.4 4 3207 763 502 20 05 2E+0  , 

0,2 9649292457 71U890E-1 

,QR- 

400 

0.4 100009 929 8695146E+0  , 

0.3 667324070 554015 3E~1 

,QR~ 

410 

0.3 7 006209 57 89 277 18 E+0  , 

0. 43ii95  08  0/659  / 6bJ8E-l 

,QR~ 

4 20 

0.3240 46825 96 848778 E+0  , 

0.4  8 8 0 9 3 26  0 5205  6 9 4 4E-1 

,QR- 

430 

0. 2 7 2710 7 3 569 4 419 7 7E+0  , 

0.5 37 22 135 057 98 28 17 E-l 

,QR- 

440 

0. 21 6896 75381 30 2257E+0  , 

0.5 7752b3 4026 8628 01E-1 

,QR- 

450 

0.1 5752133984808169E+0  , 

0.6  0835236 46 39 016 96E-1 

,QR- 

460 

0. 9 55 394 33/ 3 68081 5 0E-1  , 

0.6  291872817  3414148E-1 

,QR* 

470 

0.32028446431302813E-1  , 

0. 6396909767337 6078B-1 

/QR- 

480 

DATA  Q3 

QR- 

490 

$ 

/0 .4986319309247 4 078 E+0  , 

0. 3 509 3 050 04 / 3 50 483 E- 2 

,QR- 

500 

0.49 28 0575 57 7263 417 E+0  , 

0.813719 736 54528350E-2 

,QR- 

510 

0. 4 8 238 1127 79 3753 2 2E+0  , 

0.1269603 26 546310 30 E-l 

,QR- 

520 

0.4b745303/96B86984E+0  , 

0.171 3693 14 5551 07 17 E-l 

,QR- 

530 

0.448160 57  78 8 302606E+0  , 

0.2141794901. 113 340 B-l 

,QR- 

540 

0.4  24 6 83 80 68 66 28 49 9E+0  , 

0. 254990296311880 88E-1 

,QR* 

550 

0.3 972418979 839 71 20 E+0  , 

0. 2934 204673926 77 7 4E-1 

,QR- 

5(0 

0. 3660910593701 4484E+0  , 

0. 3291111138818 0923E-1 

,QR* 

570 

0. 3 3 152 21 33 46 510 76 0E+0  , 

0.36172897054424253 E-l 

,QR- 

580 

0.29385/H/862038116E+0  , 

0.3 9 096 94 789 353 51 53  E-l 

,QR- 

590 

0.253449954466114/0 E+0  , 

0.4 16559621134 /3378E-1 

,QR- 

600 

6 .2 106756380 6531767E+0  f 

0.4  3826046 5022019IJ6E-1 

,QR* 

610 

0.165 93 430114 10 6382E+0  , 

0.4  5 5869 393 47 8819 42E-1 

,QR- 

620 

0.1 1964 3681 126 06 854E+0  , 

0.46922199540402283E-1 

,QR- 

630 

0.72235980791398250E-1  , 

0.47819360339637430 E-l 

,QR- 

640 

0 .2 415383 284 38691 58E-1  , 

0.4  8 27 0044 25 7 3 6390 0E-1 

/QR- 

650 

$,NQ/2,4,8,12 

,12,16,24,32/, 

0R- 

6(0 

$ NS/1,3,7,15 

,25,37,53,77/ 

QR- 

(70 

DO  300  L-1,8 

UK- 

6 UN 

XF(N.EQ.NO(L) ) GO  TO  301 

UH- 

690 

CONTINUE 

Q*+- 

14  t 

WRITE (5 , 905) 

N 

QR- 

HU 

rUKMATf0  calling  PARAMETER  -',15,’  INTEGRATION  NOT  POSSIBLE'  //)QR- 

720 

RETURN 

QR- 

/J0 

CONTINUE 

QR- 

740 

NP-NS(L) 

QH- 

/5N 

NE-NP+N-i 

QR- 

/60 

AX-0.5’’ (XU+XL) 

QR- 

/?0 

BX*XU-XL 

UR* 

nut 

OVAL*  (0.  ,0.) 

QR- 

798 

CVALR=OMPLX{0. ,0. ) 

QR- 

800 

DO  350  J=NP, 

NE,2 

UR- 

810 

DX*QG(J ) *BX 

UR- 

820 

AXP»AX+DX 

QK- 

83  0 

AXM*AX-DX 

QR- 

84  0 

Ki>PsRSL+AXP* 

DEL*SINGS 

QR- 

850 

RSM-RSL+AXM* 

DEL*SINGS 

QR- 

8(0 

121 


oor.  h onnorcrrrnonor 
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CFP“CF(AXP) 

CFM“CF (AXM) 

CVAL*CVAL+QG(J+1)  MCFP+CFM) 

CVALR*CVALR+QG(J+1)  * (CFP*CMPLX  (RSP,  0 0)  KFM*CHPLX(RSM,fl.0) ) 
350  CONTINUE 

CVAJ>CVAL*BX 

CVALR=CVALR*BX 

RETURN 

END 


SUBROUTINE  ICRMUL(CM,CV  ,CVS,NMROHS, NMCOLS, NVROWS, NVCOLS) 

r w w x » ft »«  «*««****  ##r  ft*  *********  ****wwwtj*  ««  xx  xx  hr  *v  tt*w  xx  ww  ww 


PERFORMS  l CVS]  ■ l CM]  X tCVJ 

CM  * INPUT  MATRIX  (NMROWS  X NMCOLS) 

CV  - INPUT  MATRIX  (NVRCWS  X NVCOLS) 

CVS  ■ SOLUTION  MATRIX  (NMROWS  X NVCOLS) 
NMROWS  - « OF  ROWS  IN  ICM] 

NMCOLS  - « OF  COLUMNS  IN  l CM] 

NVROWS  - * OF  ROWS  IN  [CV] 

NVCOLS  - * OF  COLUMNS  IN  [CV] 


IMPLICIT  COMPLEX  (C) 

DIMENSION  CV( NVROWS, NVCOLS) ,CM(NMROWS, NMCOLS) , CVS (NMROWS, NVCOLS) 
I FINMCCLS.NE. NVROWS)  GO  TO  40 
C0-  CMPLX(0. 0,0.0) 

DO  30  J»l, NVCOLS 
DO  20  1*1, NMROWS 
CSUM-C0 

DO  13  K-l, NMCOLS 
10  CSUM*CSUM+CM(I,K) *CV(K,J) 

20  CVS(I,J) -CSUM 
30  CONTINUE 
RETURN 

40  WRITE(5, 10000)  NMCOLS, NVROWS 
00140  rUKMATC  MULTIPLICATION  NOT  POSSIBLE.*/ 

S'  NMCOLS  * ' I, 10X, 'NVROWS  - *1) 

STOP 
END 


SUBROUTINE  BESSEL  (X,N ,BJ , BV ,NX) 

IMPLICIT  REAL'S  (A-H,0-Z) 

DIMENSION  BJ (1) ,BY(1) 

REAL* 4 SNGL 
D * 1.  D-07 

CHECK  FOR  ERRORS  IN  N AND  X 

IF  (N)  710,720,720 
710  WRITE  (5,900) 

*00  FORMAT  ( 30H  N IS  NEGATIVE  ) 

GO  TO  999 

720  IF  (X)  730,730,740 


QR~ 

8/0 

QR- 

8S0 

QR- 

890 

QR~ 

900 

QR- 

910 

QR- 

920 

QR- 

930 

QR- 

940 

UR- 

950 

MM- 

10 

**MM- 

20 

MM- 

30 

MM- 

40 

MM- 

50 

MM- 

b0 

. MM- 

70 

MM- 

80 

MM- 

90 

MM- 

100 

MM- 

110 

MM- 

120 

MM- 

130 

»*MM- 

140 

Mn- 

150 

MM- 

160 

MM- 

170 

MM- 

18  0 

MM- 

190 

MM- 

200 

MM- 

210 

MM- 

220 

MM- 

230 

MM- 

240 

MM- 

250 

MM- 

2bi 

MM- 

2/0 

MM- 

28  0 

MM- 

290 

MM- 

300 

MM- 

ilB 

JY- 

10 

JY- 

20 

JY- 

J0 

JY- 

40 

J Y- 

50 

JY- 

60 

JY« 

70 

JY- 

80 

JY- 

90 

JY- 

100 

JY- 

119 

JY- 

120 

JY- 

130 

122 


-fel&FAQflt  IS  BIST  QUALITY  PRACTICABLB 
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730 

WRITE  (5,910) 

JY- 

140 

010 

FORMAT  (30K  X IS  NEGATIVE  OR  ZERO 

) 

JY- 

150 

GO  TO  999 

JY- 

160 

740 

IF  (NX  .EQ.  1 ) GO  TO  30 

JY- 

1/0 

IF  (NX  .EQ.  2)  GO  TO  600 

JY- 

180 

IF  (NX  ,.EQ.  3)  GO  TO  30 

JY- 

190 

WRITE  (5,920) 

JY- 

200 

020 

FORMAT  (30H  ERROR  IN  PARAMETER  NX 

) 

JY- 

210 

000 

WRITE  (5,930)  X, N ,NX 

OY- 

220 

030 

FORMAT  (30H  FARAMKTER  VALUES  X,N,NX 

,G15.5 ,215  ) 

JY- 

230 

STOP 

JY- 

240 

C 

JY- 

250 

C CALCULATION  OF  J BESSEL  FUNCTION 

JY- 

260 

C 

JY- 

2/0 

30 

L"N 

JY- 

280 

31 

IF(X-15. ) 32,32,34 

JY- 

290 

32 

NTEST“20.+10.*X-X**  2.  / 3. 

JY- 

300 

GO  TO  36 

. JY- 

310 

34 

NTEST-90.+X/2. 

JY- 

320 

36 

IF(L-NTEST) 40,38,38 

JY- 

3 30 

38 

WRITE  t 5,940) 

JY- 

340 

040 

FORMAT  (3  0H  RANGE  OF  X VIOLATED 

) 

JY- 

350 

GO  TO  999 

J Y“ 

360 

40 

N1  - L + 1 

JY- 

370 

BJ(N1) -0.0 

JY- 

380 

BPREV-.0 

JY- 

390 

C 

JY- 

400 

C 

COMPUTE  STARTING  VALUE  OF  M 

JY- 

410 

C 

JY- 

420 

IF(X-5. )50,60,60 

JY- 

4 30 

50 

MA-X+6. 

J X- 

440 

GO  TO70 

J X- 

450 

60 

MA-1 . 4 "X+60./X 

JY- 

460 

70 

Mb  * L + i FIX  ( SNGL  (X)  )/4  +2 

JY- 

4/0 

M ZERO- MAX 0 (MA,MB) 

JY- 

480 

C 

JY- 

490 

C 

SET  UPE\9L\\\\TOo6MM 

JY- 

500 

C 

JY- 

510 

MMAX-NTEST 

JY- 

5 20 

100 

DO  190  M«MZERO,MMAX,3 

JY- 

530 

C 

JY- 

540 

C 

SET  F(M) , F(M-l) 

JY- 

550 

C 

JY- 

560 

FM1-1.0E-28 

J X- 

5 /0 

FM-.0 

JY- 

580 

ALPHA-. 0 

J Y- 

5 90 

IF(M- (M/2) *2)120,110,120 

UY- 

600 

110 

JT--1 

JY- 

bl0 

GO  TO  130 

JY- 

620 

120 

JT-1 

JY- 

630 

130 

M2-M-2 

JY- 

6 40 

DO  160  K-1.M2 

jy- 

650 

MK-M-K 

JY- 

660 

BMK  * 2. "DFLOAT(MK) *FM1/X-FM 

JY- 

670 

FM-FM1 

JY- 

680 

FMl-BMK 

JY- 

6 90 

IF(MK-L-i) 150,140,150 

JY- 

70  0 

140 

BJ(L4-1)  “BMK 

JY- 

710 

150 

JT--JT 

J Y- 

/20 

S-1+-JT 

J X- 

/3  0 

160 

ALPHA“ALPHA+BMK*S 

JY- 

740 

12° 
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BMK*2.*FMl/X-FM 
IF(L)  180, 170, 3.80 
170  BJ(L+1)=BMX 
180  ALPHA=ALPHA+BMK 

B J ( L+ 1 ) *BJ ( L+ 1)  /ALPHA 
ERR*BJ (L+l) -BPREV 

IF (DABS (ERR) -DABS(D*BJ(L+1) ) ) 200/200,185 
185  BPREV*BJ (L+l) 

190  ERR“ERR/BJ(L+1) *100. 

SSU  rOKMAT(59H  KtSUUlKED  ACCURACY  IN  BJN  NOT  OBTAINED, 
US  ,G20.6) 

200  IF(L-N)210, 220,220 

220  L**L-1 

IF(L  .LT.  0)GO  TO  240 
GO  TO  31 

210  IP(L  .EQ.  0>UO  TO  240 
.00  230  I-2,N 
1>K-I 

230  1JJ  (L+i)  "*2.*DFLOAT(L+l)  *BJ ( L+2)  A-B J(L+3) 

240  IF  (NX  .EQ.  3)  GO  TO  600 
RETURN 
C 

C CALCULATION  OF  Y BESSEL  FUNCTION 
C 

600  PI  ■ 3.141592653 


BRANCH  IF  X LESS  THAN  OR  EQUAL  4 


JY- 

75  f 

JY- 

7 66 

JY- 

770 

JY- 

780 

JY- 

790 

JY- 

800 

JY- 

810 

JY- 

820 

JY- 

83  0 

JY- 

840 

PER  CENT  DIFF.  JY- 

850 

JY- 

8fc0 

JY- 

8/0 

JY- 

880 

JY- 

390 

J Y- 

900 

JY- 

910 

JY- 

920 

JY- 

930 

JY- 

940 

JY- 

950 

JY- 

960 

JY- 

9/0 

JY- 

980 

JY- 

990 

JY- 

1000 

JY- 

1010 

JY- 

1020 

JY- 

1030 

IF  ( X— 4 . ) 640,640,630 

e COMPUTE  Y0  AND  Y1  FOR  X GREATER  THAN  4 

C 

630  T-4./X 

P0-. 3989422793 
Q0*-.0 124669441 
P1-. 3989422819 
Ql*. 0374008364 
A«T*T 


JY-  1040 
JY-  1 050 
JY-  1068 
JY-  1070 
JY-  1080 
JY-  1090 
JY-  1100 
JY-  1118 
JY-  1120 
JY-  1130 
JY-  1140 


B-A 

JY-  1150 

P0»P0-.0017530620*A 

JY-  1160 

Q0-Q0+. 0004564 324*A 

JY-  1176 

P1-P1+.029218256*A 

JY-  1180 

Q1«Q1-.00063904*A 

JY-  1198 

A-A*A 

JY-  1200 

P0*P0+v 00017343* A 

JY-  1210 

Q0*‘Q0-.0C00869791*A 

JY-  1220 

Pl*Pl-. 000223 20 3*A 

JY-  1230 

Q1*Q1+. 8001064741* A 

JY-  124.) 

A«A*B 

JY-  1250 

P0*P 0-. 00004876 13*A 

JY-  1260 

Q0*Q0+. 000 0342468* A 

JY-  127  0 

Pl*Pl+.0 800580759 *A 

JY-  1280 

Q1»Q1-.0000398708*A 

JY-  1290 

A»A*B 

JY-  1300 

P0*P0+.0000173565*A 

JY-  1310 

Q0»Q8-.0 0J0142078 *A 

JY-  1320 

PI  “PI-. 000020092 *A 

JY-  1330 

01*01+. 000 015 2 2* A 

jy-  1340 

A=A*B 

P0»P0~,8000037043*A 

JY-  1350 

124 
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Q0-Q0+. 000003 2J12*A 
Pl3,Pl+.  0 0000  4 2414  *A 
Qi*=Ql-.0000036594*A 
A * DSQRT(2.*P1) 

B«4.*A 

P0*A*P0 

Q0-B*Q0/X 

P1»A*P1 

Ql«B*Ql/X 

A»X-P 1/4. 

B * DSQRT  { 2./  (PI*X) ) 

Y0  * B* (P0*DSIN(A) +Q0*DCOS(A) ) 

Yl  - 3* (-Pl*DCOS(A) +Q1*DSIN(A) ) 

GO  TO  690 
C 

C COMPOTE  Y0  AND  Yl  FOR  X LESS  THAN  OR  EQUAL  TO  4 

C 

640  XX-X/2. 

X2*XX*XX 

T-DLOG(XX) +.5772156649 
SUM-0. 

TERM-T 

Y0-T 

DU  670  L - 1,15 
IF(L-l)  650,660,650 
650  SUM^SUM+l./ OFLOAT(L-l) 

660  FL-L 

TS-T-SUM 

TERM- (TERM* (-X2) /FL*«2) * (l.-l./ (FL*TS) ) 

670  Y0-Y0+TERK 

TERM  - XX* (T-  ,5 } 

SUM-0. 

Yi*TERM 

DO  600  L - 2,16 
SUM-SUM+l./DFLOAT(L-l) 

FL-L 

FL1-FL-1. 

TS-T-SUM 

TERM- (TERM* (-X2) /(FL1*FL) ) * ( (TS- ,5/FL) / (TS+.5/FL1) ) 
600  Y1-Y1+TERM 
PI2-2./P1 
Y0-P 12" Y 0 
Y1--F12/ X+P12"Y1 
C 

C CHECK  IF  ONLY  Y0  OR  Yl  IS  DESIRED 
C 

600  XF(N-l)  500,500,530 
C 

v.  Kfci'UKN  E.TtlEK  Y0  OK  Yl  AS  KEUUIKEU 

v. 

500  1 F ( N } 510,520,510 

510  BY ( 2 ) -Y  1 
520  BY ( 1 ) -Y  0 

RETURN 

C 


JY-  1360 
JY-  1370 
JY-  1380 
JY-  1390 
JY-  1400 
JY-  1410 
JY-  1420 
JY-  1430 
JY-  1440 
JY-  1450 
JY-  1460 
JY-  147C 
JY-  1480 
JY-  .1490 
JY-  1500 
JY-  1510 
JY-  1520 
JY-  1530 
JY-  1540 
JY-  1550 
JY-  1560 
JY-  1570 
JY-  1580 
JY-  13*0 
JY-  1600 
JY-  1610 
JY-  1620 
JY-  1630 
JY-  1640 
JY-  1650 
JY-  1660 
JY*  16/0 
JY-  1680 
JY-  1650 
JY-  1700 
JY-  1/li 
JY-  1720 
JY-  1/30 
JY-  1/40 
JY-  17  50 
JY-  1760 
JY-  1 / /0 
JY-  1/80 
JY-  1/90 
JY-  18  00 
JY-  1810 
JY-  1820 
Jl-  163  0 
JY-  104  0 
JY-  10  50 
JY-  1060 
JY-  10/0 
JY-  1880 
JY-  1890 

Jl-  190  6 


t-BKrUKM  KbCUKXENCE  OPERATIONS  TO  FIND  YN'.X)  JY-  1910 

C J X-  1 92  0 

530  BYU;-xe  JY-  1930 

BY ( 2 j -Y 1 JY-  1940 

DO  545  K-2,N  JY-  1950 

T=*DFLOAT(2*  (K-l)  )/X  J Y-  i960 


12 


Sf  15  BBST  men  casl* 

JWaiSHBD  ?0  DDQ 


BY  (K+l)  «T*BY(H)  -BY(K-l) 

JX-  19/0 

c 

JX*  1 9019 

V 

BESSEL  X FUNCTION  HAS  EXCEEDEU  119**35 

JY-  1990 

c 

JX-  2000 

IP  ( DABS(BY(K+1) )•  1.0D35) 

545,545,541 

JY-  2 BIB 

54] 

WHITE  (5/980) 

JX-  2020 

980 

tUKMAT  (3t)H  X BESSEL  FCN 

EX  CEEDED  10**35  ) 

JY-  2030 

GO  TO  999 

JX-  2B4B 

545 

CONTINUE 

JY-  2050 

KETUKN 

0 1-  2 000 

END 

oi-  2»/» 

FUNCTION  ELIClK(AMl) 

CE- 

IB 

C** 

20 

c 

CE- 

3B 

c 

COMPLETE  ELLIPTIC  INTEGRAL  OP  THE 

FIRST  KIND  K(M),  . CE- 

40 

c 

WHERE  AMl-l-M 

CE- 

50 

c 

CE- 

60 

c 

REFERENCES 

CE- 

70 

c 

CE- 

80 

c 

ABRAMOWITZ  AND  STEGUN,  EQ.  17 

.3.34  CE- 

90 

c 

CE- 

100 

c 

MAGNITUDE  (ERROR)  .LE.  2.0E-8 

CE- 

110 

c 

CE- 

120 

C” 

DATA  A0 , A1 , A2  , A3  , A4  , B0 , 81 , B2  , B3 , B4/ 

CE- 

130 

140 

5 1.3U62943bil2,  .09666344259, 

.03590092383,  .03742563713,  CE- 

150 

$ .01451196212,  .5, 

.12498593597,  .06880248576,  CE- 

160 

$ .03328355346,  .00441787012/ 

CE- 

l/B 

A-A0+A1*AM1 

CE- 

100 

B-B0+B1*AM1 

CE- 

190 

IP (AMI  .LT.  l.E-18)  GO  TO  10 

CE- 

200 

AMI 2“AMl* AMI 

CE- 

210 

A-A+A2-AM12 

CE- 

220 

B-B+B2*AM12 

CE- 

230 

IP(AM1  .LT.  l.E-12)  GO  TO  10 

CE- 

240 

AM13-AM1  2*  AMI 

CE- 

250 

A»A+A3"AM13 

CE- 

2 60 

B-B+BJ*AM13 

CE- 

270 

IP (AMI  .LT.  i.E-9)  GO  TO  10 

CE- 

2 80 

AM14«AM13*AM1 

CE- 

230 

A"Ai-A4»AHl4 

CE- 

300 

B«B+B4*AM14 

CE- 

310 

IB  CONTINUE 

CE- 

3 20 

ELIC1K=*A-B*ALQG(AM1) 

CE- 

3 30 

RETURN 

CE- 

340 

END 

Cis- 

350 

subroutine  TRPADP(CF,XL0W,XH1GH,ER,MAXP,IER,CANS) 

AT- 

10 

c*w 

*A1— 

20 

C 

AT- 

30 

c 

INTEGRATION  OF  A FUNCTION  BY  ADAPTIVE  TRAPEZOIDAL  RULE 

AT- 

40 

C 

AT*- 

50 

c 

AT- 

60 

c 

CE  = EXTERNALLY  SUPPLIED  COMPLEX  FUNCTION  TO  BE  INTEGRATED 

AT- 

70 

C 

AT- 

80 

XLOW  = LOWER  LIMIT  OF  INTEGRATION 

AT- 

90 

126 
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C A'l- 
C XHIGH  » ({UPPER  LIMIT  OF  INTEGRATION  MINUS  LOWER  LIMIT)  DIVIDED  AT- 
C BY  TWO)  + LOWER  LIMIT  AT- 
C AT” 
C ER  RELATIVE  ERROR  CRITERION  FOR  CONVERGENCE  CHECK  Al- 
C AT- 
C MAXP  ■ MAXIMUM  NUMBER  OF  POINTS  TO  BE  USED  AT- 
C AT- 
C IGR  **  ERROR  CODGs  AT“ 
C * 1,  IF  INTEGRATION  DID  NOT  CONVERSE  WITHIN  (MAXP)  AT- 
C POINTS  AT- 
C “ 0/  OTHERWISE  Al- 
C A*iw 
C CANS  « RESULT  OF  THE  INTEGRATION  AT- 
<;  AT- 
C AT” 
C CF  MUST  BE  LISTED  IN  AN  EXTERNAL  STATEMENT  IN  THE  AT- 
C CALLING  PROGRAM.  AT- 


C 


AT- 


«- 

c 

c 

c 

c 


THE  RESULT  (CANS)  IS  INTEGRAL (CF)  BETWEEN  THE  LIMITS  (XLOW)  AND  AT- 
(XLOW+ (XHIGH-XLOW) *2) . SYMMETRY  OF  THE  FUNCTION  (CF)  IS  ASSUMED  AT- 


ABOUT  THE  POINT  (XHIGH) . 

AT- 

AT” 

IMPLICIT  CONPLEX(C) 

AT- 

COMMON/NPOINT/NP 

AT- 

IER-0 

AT- 

DELO« XHIGH-XLOW 

AT- 

NPOLD-2 

AT- 

C2«CMPLX(2. 0,0.0) 

AT- 

IP(MAXP  .EQ.  0)  MAXP-10000 

AT- 

IF(ER  .EQ.  0)  EI^l.E-4 

AT- 

CDELOCMPLX(  DELO, 0,0 ) 

AT- 

COLD- (CF(XLOW) +CF (XHIGH) ) /C2 

AT- 

C-COLD 

AT- 

NP-N  POLD+N  ROLD-1 

AT- 

DEL-DELO/2.0 

AT*- 

CDEL-CMPLX(D8L,0.?) 

AT- 

K-NP-NPOLD 

AT- 

K-K+K 

AT- 

DO  20  J-2,K,2 

AT- 

X“XLOW+DEL*FLOAT(J-l) 

AT- 

C-C+CFU) 

AT- 

CO-COLD*CDELO 

AT- 

CN-C*CD3L 

AT- 

CND-CN 

AT- 

ABSND-CABS (CND) 

AT- 

IFfABSND  ,LT.  l.E-20)  ABSND«i.  B-20 

AT- 

DIFF-CABS(CN-CO) 

AT- 

ERR-DIFF/ABSND 

AT- 

IF (ERR  .LT.  ER)  GO  TO  40 

AT- 

IF (NP  .GE.  MAXP)  GO  TO  30 

AT- 

DELO-DEL 

AT- 

CDELO-CDEL 

AT- 

NPOLD-NP 

AT- 

COLD-C 

AT- 

GO  TO  10 

AT- 

IER-1 

AT- 

CANS«CN*C2 

Ai- 

RETURN 

AT- 

END 

AT- 

IBB 

110 

120 

1 10 
140 
150 
160 
1/0 
IBB 
190 
200 
210 
220 
210 

2 40 

2 50 
260 
270 

250 

290 

300 

110 

120 

110 

3 40 
150 
360 

3 70 
1H0 
19  0 
400 
410 

4 20 
430 
4 40 
450 
4 60 
4 70 
4 00 

4 90 
500 
510 
520 
530 

5 40 
550 
560 
570 
5B0 
590 
600 
610 
620 
610 
64  0 
650 
660 
6/0 
5B0 
690 
/00 
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